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1.0  Introduction  and  Summary 

The  objective  of  this  project  is  to  optimize  surface  wave  measurements,  particularly  at  regional 
and  local  distances,  and  over  the  high  frequency  period  band  of  8-15  seconds.  There  are  three 
main  parts  to  this  study:  1)  We  collect  and  then  invert  Eurasian  surface  wave  amplitude  data  for 
attenuation  and  Q  structure;  2)  we  develop  amplitude  corrections  using  the  Bom  approximation 
and  then  compare  them  with  observations  and  with  the  results  of  finite  difference  calculations;  3) 
we  develop  a  time  domain  path  corrected  surface  wave  magnitude  and  compare  it  with  other 
surface  wave  magnitudes.  In  addition,  we  analyze  the  surface  waves  from  the  North  Korean 
nuclear  test  of  October  9,  2006.  The  path  corrected  surface  wave  magnitude  is  similar  to  the 
Butterworth  filtered  magnitude  developed  by  Russell,  but  regionalized  to  take  into  account 
variations  in  earth  stmcture  and  attenuation.  We  conclude  that  the  Bom  approximation  is  of 
marginal  value  for  attenuation  studies  because  it  is  a  small  perturbation  approximation  and  the 
observed  amplitude  variations  are  too  large  and  complicated  to  be  predicted  by  the  Born 
approximation.  We  find  that  the  average  surface  wave  Q  of  Eurasia  is  lower  than  our  previous 
background  model,  particularly  along  a  band  mnning  through  the  Middle  East.  We  develop 
attenuation  maps  that  predict  attenuation  along  any  path  in  Eurasia  for  frequencies  below  0.15 
Hz  for  distribution  to  AFRL,  the  DOE  Knowledge  Base  and  other  researchers. 

1.1  Attenuation 

We  inverted  surface  wave  amplitude  data  for  attenuation  and  corrections  to  source  moment  for 
data  from  about  300  Eurasian  earthquakes.  We  used  two  data  sets:  one  our  own  measurements 
and  one  set  from  Anatoli  Levshin  at  the  University  of  Colorado.  The  data  sets  are  fairly 
consistent  and  both  indicate  higher  attenuation  than  our  initial  background  model.  There  is  a 
large  amount  of  scatter  in  both  data  sets,  likely  the  result  of  stractural  variations  and  interference. 
Nevertheless,  there  is  enough  redundancy  in  the  data  that  we  were  able  to  perform  Q  inversions 
for  the  Eurasian  continent.  Figure  1.1  shows  a  map  of  the  attenuation  coefficients  determined  by 
the  inversion  at  10  seconds  period,  and  the  data  fit  averaged  over  all  paths. 
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Figure  1.1.  Attenuation  coefficients  at  10  seconds  period  in  Eurasia  derived  by  inversion  of  surface  wave 
amplitude  measurements  (left),  average  data  and  average  data  fit  (right).  Blue  curve  shows  the 
inversion  results,  the  red  curve  the  average  of  the  data,  and  the  green  line  the  starting  model, 
averaged  over  all  paths  longer  than  5000  km.  Dashed  lines  show  one  standard  deviation  above  and 
below  each  curve. 
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1.2  Born  approximation 

We  modeled  the  effect  of  heterogeneous  structure  using  the  Bom  approximation.  We  performed 
both  data  analysis  and  large  3-dimensional  finite  difference  calculations  to  assess  the 
performance  of  the  Bom  approximation.  Our  goal  was  to  use  the  Born  approximation  to  correct 
for  scattering  and  diffraction  caused  by  heterogeneous  stracture  prior  to  performing  Q 
inversions.  However,  the  stmctural  complexity  appears  to  exceed  the  limits  of  the  Bom 
approximation,  and  application  of  the  Bom  corrections  do  not  improve  inversion  results  at  these 
high  frequencies.  Figure  1.2  shows  a  comparison  of  the  Bom  approximation  with  the  results  of  a 
3D  finite  difference  calculation  at  10  seconds  period.  Both  model  an  explosion  near  the  right  side 
of  the  grid  in  a  background  medium  typical  of  Eurasian  shield  regions,  with  an  inclusion  of  low 
velocity  material  modeled  after  the  Tarim  Basin.  Although  there  is  similarity  between  the  two 
calculations,  the  Bom  approximation  does  not  reproduce  the  fine  stmcture  of  the  numerical 
calculations,  and  the  disagreement  is  large  enough  that  the  Bom  corrections  are  inadequate  for 
the  purposes  of  correcting  surface  wave  amplitudes  at  individual  receiver  points. 
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Figure  1.2.  Comparison  of  Born  approximation  (left)  with  finite  difference  calculation  (right)  of  amplitude 
perturbations  at  10  seconds.  The  rectangular  inclusion  is  modeled  after  the  Tarim  Basin  structure, 
and  the  external  structure  after  a  Eurasian  shield  earth  structure.  The  source  is  on  the  horizontal 
axis  at  the  right  edge  of  the  plot.  The  amplitude  is  increased  in  a  band  above  and  to  the  left  of  the 
inclusion  in  both  cases,  and  decreased  above  that.  However,  there  are  some  interference  effects  in  the 
finite  difference  calculation  that  are  not  reproduced  in  the  Born  calculation. 


1.3  Path  corrected  magnitude 

Russell  (2006)  proposed  a  new  type  of  surface  wave  magnitude  Ms(b)  that  uses  a  Butterworth 
filter  to  measure  a  time  domain  amplitude  in  a  narrow  band  around  any  desired  frequency,  and 
then  applies  a  correction  for  the  frequency  dependence  of  an  explosion  source  function.  The 
main  purpose  of  Ms(b)  is  to  allow  surface  waves  to  be  measured  at  regional  distances  at  higher 
frequencies  than  traditional  20  second  Mj.  The  magnitude  is  defined  by 


^,(6)  =log(4)  +  ^log(sinA)  +  0.0031 


r2o 

\T  ) 


A -0.66  log  —  -log  (X)- 0.43 
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where  Ab  is  the  filtered  amplitude,  T  is  the  measured  period,  and  fc  is  the  Butterworth  filter 
width. 

A  regionalizeable  path  corrected  time  domain  magnitude  can  be  derived  by  combining  the  path 
corrected  spectral  magnitude  (Stevens  and  McLaughlin,  2001)  with  Ms(b),  using  the  source  and 
path  corrections  from  earth  models  to  replace  the  empirical  corrections.  We  define  the  path 
corrected  time  domain  magnitude  Ms(bp)  as: 

=  >og  (  A  )  +  ^  log  (sin  A)  +  log  e  -  log  (5, )  -  log  (^^ )  -  log  (/J  +  (1 .2) 

Where  Si  and  S2  are  source  and  path  excitation  functions,  jp  is  the  attenuation  function  between 
the  source  and  receiver,  and  Cbp  is  a  constant  chosen  to  make  Ms{bp)  consistent  with  historical 
magnitude,  which  is  shown  to  be  Cbp=-17.96. 

The  Russell  Butterworth  filtered  magnitude  and  the  path  corrected  time  and  frequency  domain 
magnitudes  have  a  similar  purpose,  specifically  to  allow  surface  waves  to  be  measured  at 
regional  and  local  distances  and  at  higher  frequencies  than  20  seconds.  These  magnitudes  give 
similar  results  when  applied  to  data  and  are  consistent  in  value  with  traditional  20  second 
magnitudes.  The  path  corrected  magnitudes  have  the  advantage  that  they  can  be  regionalized  to 
take  into  account  differences  in  earth  structure  and  attenuation,  while  the  Butterworth  filtered 
magnitude  uses  a  good  representative  average  value  for  these  quantities.  An  issue  with  all  of  the 
magnitudes  is  how  to  determine  which  frequency(ies)  to  measure.  The  path  corrected  spectral 
magnitude,  for  example,  performs  a  robust  average  over  all  frequencies,  while  the  Bonner  et  al. 
implementation  of  the  Butterworth  filtered  magnitude  uses  the  maximum  value.  Using  the 
maximum  value  may  give  better  discrimination  but  at  the  cost  of  more  variability  in  the 
magnitude. 

1.4  Korea 

The  North  Korean  nuclear  test  of  October  9,  2006  generated  larger  than  expected  surface 
waves,  nearly  a  magnitude  unit  higher  than  would  be  expected  based  on  larger  events 
assuming  a  Ms:log(yield)  slope  of  one.  Some  of  this  difference  can  be  explained  by  the 
absence  of  tectonic  release  and  high  velocity  source  medium  for  this  event.  However,  the 
results  also  suggest  that  the  Ms:log(yield)  slope  may  be  less  than  one,  which  has  implications 
for  discrimination  of  small  events  as  well  as  yield  estimation.  Figure  1.3  shows  path 
corrected  magnitudes  for  7  stations  that  recorded  the  North  Korean  explosion.  Although 
there  is  some  variation  due  to  earth  structure,  they  clearly  show  an  Ms  of  about  2.9.  The 
right  side  of  Figure  1.3  shows  a  fit  to  the  Mjryield  data  for  all  other  explosions  with  known, 
unclassified  yield  and  Ms,  showing  that  the  predicted  Ms  for  this  event  is  approximately 
2.0.  However,  this  is  based  on  a  slope  of  1  in  the  Ms:log(yield)  curve,  which  is  expected 
from  cube  root  scaling,  but  may  not  be  right.  The  other  curves  show  Mueller-Murphy 
source  predictions  for  standard  containment  depth  and  for  an  overburied  explosion  (believed 
to  be  the  case  here),  and  those  curves  are  much  more  consistent  with  the  North  Korean 
data  point. 
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M  vs.  Yield 


Figure  1.3.  Path  corrected  surface  wave  magnitudes  for  the  North  Korean  nuclear  test,  and  fit  to  three 
MsiYield  curves  for  an  estimated  yield  of  1  kiloton,  together  with  historical  MsiYield  data. 

1.5  Recommendations 

Additional  research  is  needed  in  the  following  areas: 

1.  The  high  Ms  of  the  North  Korean  earthquake  suggests  that  the  Msilog  yield  slope  may  be 
less  than  one,  which  has  implications  for  discrimination  of  small  events  as  well  as  yield 
estimation.  More  research  is  needed  to  analyze  surface  waves  from  small  events  of 
known  yield  in  high  velocity  media  to  see  if  the  North  Korean  event  is  an  anomaly  or  the 
norm. 

2.  The  Butterworth  fdtered  magnitudes  and  spectral  magnitudes  are  all  sensitive  to  how  the 
frequency  for  measurement  is  picked.  As  discussed  above,  using  the  maximum  value  may 
give  better  discrimination  but  at  the  cost  of  more  variability  in  the  magnitude,  and  that 
advantage  may  go  away  for  small  events.  More  research  is  needed  to  determine  the 
optimum  frequency  strategy. 

3.  A  method  for  correcting  surface  wave  amplitudes  for  path  structure  is  needed.  As 
discussed  above,  the  Bom  approximation  is  not  adequate  for  stmctures  as  complex  as 
Eurasia  at  frequencies  of  8-15  seconds.  Other  approximation  methods  have  similar 
problems,  and  a  new  procedure  may  be  required. 

4.  The  robustness  of  the  Q  models  for  Eurasia  depend  on  the  amount  of  data  that  goes  into 
the  inversion.  That  is  particularly  trae  in  the  absence  of  an  adequate  method  for 
correcting  for  earth  stmcture.  So  we  recommend  that  a  much  larger  data  set  be  processed 
and  included  in  the  Q  inversion. 
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2.0  Path  Corrected  and  Frequency  Dependent  Surface  Wave  Magnitudes 

Surface  wave  magnitudes  play  an  important  role  in  earthquake/explosion  discrimination.  Three 
main  problems  exist  with  traditional  surface  wave  magnitudes:  1)  surfaee  wave  dispersion  causes 
amplitude  variations  unrelated  to  the  source;  2)  it  is  not  possible  to  reliably  measure  a  traditional 
20  second  surface  wave  magnitude  at  local  and  regional  distances  because  the  surface  wave  is 
not  dispersed  enough;  and  3)  differences  in  earth  structure  and  attenuation  cause  variations  in 
surface  wave  amplitudes  that  are  unrelated  to  the  source.  Several  surface  wave  magnitude 
measurements  have  been  proposed  to  address  these  limitations,  and  we  propose  some  further 
improvements  in  the  form  of  a  regionalized  path  corrected  surface  wave  magnitude.  In  the 
following  section  we  compare  these  magnitude  types  in  detail,  and  then  apply  them  to  the  North 
Korean  nuclear  test.  We  also  use  this  as  an  example  of  separating  the  surface  waves  into  their 
component  source,  receiver  and  path  parts,  and  examine  the  residual  to  show  the  effect  of 
structural  variations. 


2.1  Comparison  of  Butterworth  filtered  and  regionalized  surface  wave  magnitudes 

Russell  (2006)  proposed  a  new  type  of  surface  wave  magnitude  Ms(b)  that  uses  a  Butterworth 
fdter  to  measure  a  time  domain  amplitude  in  a  narrow  band  around  any  desired  frequeney,  and 
then  applies  a  correction  for  the  frequency  dependence  of  an  explosion  source  function.  The 
main  purpose  of  Ms(b)  is  to  allow  surface  waves  to  be  measured  at  regional  distances  at  shorter 
periods  than  traditional  20  second  Mj.  The  magnitude  is  defined  by 


=  log  (4 )  +  ^  log  (sin  A)  +  0.003 1  ^ 


V  J  y 


a'-*  ( 20^ 

A -0.66  log  —  -log  (/J- 0.43 


(2.1) 


where  Ab  is  the  filtered  amplitude,  T  is  the  measured  period,  and  f  is  the  Butterworth  filter 
width.  This  magnitude  also  requires  that  the  frequency  band  be  less  than  a  minimum  value 

Q 

defined  by  f  <  — ^ .  Russell  (2006)  finds  Gmin=0.6  for  continental  structures  at  periods 
rVA 

between  8  and  40  seconds,  with  smaller  values  required  for  deep  sediment  structures.  Russell 


(2006)  also  shows  that  4  = 


A  where  G  is  a  constant  which  for  typical  continental  paths 


is  approximately  0.93,  and  A  is  the  equivalent  time  domain  amplitude.  Note  that  if  Gmin  is  fixed, 
then  the  filter  correction  corresponds  to  a  distance  correction  for  a  normally  dispersed  (non-Airy 

phase)  surface  wave  of^log  A  . 

Stevens  and  McLaughlin  (2001)  defined  a  path  correeted  spectral  magnitude,  which  similarly  is 
intended  to  allow  surfaee  waves  to  be  measured  at  all  distances  and  frequencies,  and  in  addition 
is  regionalizeable  since  it  is  derived  from  earth  models.  The  path  correeted  spectral  magnitude, 
logMo',  is  calculated  by  dividing  the  observed  surface  wave  spectrum  by  the  Green’s  function  for 
an  explosion  of  unit  moment  and  taking  the  logarithm  of  this  ratio,  averaged  over  any  desired 
frequeney  band.  The  path  corrected  spectral  magnitude  is  defined  as  the  logarithm  of 


Mo  = 


U^(co,r,0)i 


^4(6i,4)4(®)exp[-y  r«)r]^ 


yja^sm(r/aj 


(2.2) 
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where  Uz  is  the  observed  vertical  component  surface  wave  spectrum,  St  depends  on  the  source 
region  elastic  structure  and  the  explosion  source  depth,  S2  depends  on  the  receiver  region  elastic 
structure,  and  jp  is  the  attenuation  coefficient  that  depends  on  the  attenuation  integrated  over  the 
path  between  the  source  and  receiver.  All  of  the  functions  in  equation  2.2  can  be  derived  from 
plane-layered  earth  models  (see  Stevens  and  McLaughlin,  2001),  and  allow  the  measurement  to 
be  regionalized  to  account  for  differences  in  earth  structure  at  the  source  and  receiver,  and  due  to 
attenuation  along  the  path.  Since  Mo'  is  a  physical  quantity,  equation  2.2  is  assumed  to  be  in  SI 
units  and  logMo'  is  in  log(nt-m),  however  in  the  following  comparison  we  express  Uz  in  nm-s  for 
consistency  with  the  other  amplitude  measurements.  This  adds  a  constant  value  of  -9  to  the 
normalization  constant  for  logMo'. 

Since  equations  2.1  and  2.2  are  both  intended  to  flatten  the  surface  wave  spectrum,  in  principle 
they  can  be  measured  over  any  desired  frequency  band.  In  practice,  the  path  corrected  spectral 
magnitude  (equation  2.2)  has  been  calculated  by  averaging  over  a  frequency  band  designed  to 
avoid  noise  contamination,  with  implementation  of  an  outlier  rejection  scheme  to  minimize  bias 
from  spectral  dips  and  noise  (Stevens  et  al.,  2005).  The  implementation  of  the  Butterworth 
filtered  magnitude  (equation  2.1)  by  Bonner  et  al.  (2006)  has  instead  used  the  maximum  value 
over  a  period  band  of  8-25  seconds,  defined  as  Ms(VMAX),  with  analyst  rejection  of  outliers. 

A  path  corrected  time  domain  magnitude  can  be  derived  by  combining  the  path  corrected 
spectral  magnitude  with  Ms(b),  using  the  source  and  path  corrections  from  earth  models  to  replace 
the  empirical  corrections.  We  define  the  path  corrected  time  domain  magnitude  Ms(bp)  as: 

=  log  (4  )  + 1  log  (sin  A)  +  A  log  e  -  log  (  Aj )  -  log  (  Aj )  -  log  (/^ )  +  (2.3) 


where  Cbp  is  a  constant  chosen  to  make  Ms(bp)  consistent  with  historical  magnitudes.  By  defining 
Ms(bp)  to  be  equal  to  Rezapour  and  Pearce  (1998)  Ms  at  50“,  and  simultaneously  using  the 
Rezapour  and  Pearce  attenuation  rate,  and  using  the  mean  20  second  value  of  Si  and  S2  for 
Central  Asian  continental  structures  (log(Si)-i-log(S2)=- 17.41),  we  find  Cbp=-17.96.  We  can  also 
define  a  spectral  magnitude  directly  from  equation  2.2,  using  the  relation  (again  from  Russell, 
Att 


2006) 


A=- 


fJJ^.  This  gives  logM^ =11.74,  which  is  identical  to  the  mean 


difference  between  logMo'  and  Rezapour  and  Pearce  Ms  found  through  measurement  of  a  large 
data  set  by  Stevens  and  McLaughlin  (2001).  We  can  therefore  define  an  equivalent  spectral  Mg, 
which  we  define  as  Ms(sp)=logMo'- 11.74,  which  adding  the  logMo'  normalization  constant 
V2log(ae)-9  gives  an  Ms(sp)  normalization  constant  of  -17.34.  Table  2.1  shows  a  comparison  of  the 
terms  in  each  of  these  magnitudes,  and  in  the  Rezapour  and  Pearce  Mg. 
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Table  2.1.  Comparison  of  time  domain  and  spectral  magnitude  measurement  and  correction  terms 


Magnitude 

Type 

Amplitude 

Measure 

Source 

Receiver 

Geometric 

Spreading 

Attenuation 

Dispersion 

Filter 

Norm 

M, 

log(A/T) 

—  log  (sin  a) 
2 

0.0046A 

1 

-log  A 

3 

2.37 

-l^s(b) 

log(Ab) 

f  20 

-0.66  log  — 

\T  J 

—  log(sin  a) 
2 

f20Y* 

.0031  —  A 

VT  J 

-0.43 

logMo' 

log(Ud 

-log(S,) 

-l0g(S2) 

—  log(sin  a) 
2 

y^Aloge 

—  log  a  —9 

2 

-l^s(sp) 

log(Ud 

-log(S,) 

-l0g(S2) 

—  log(sin  a) 
2 

y^Aloge 

-17.34 

^^s(bp) 

log(Ab) 

-log(Si) 

-l0g(S2) 

—  log  (sin  a) 
2 

y^Aloge 

-17.96 

In  this  table  A  is  the  traditional  time  domain  20  second  amplitude  in  nm,  Ab  is  the  Butterworth 
fdtered  magnitude  (using  a  3  pole  two  pass  phaseless  fdter)  in  nm,  and  Uz  is  the  Fourier  spectral 
amplitude  in  nm-s.  Figure  2.1  (left)  shows  a  comparison  of  Russell’s  approximation  to  the 
explosion  excitation  function  with  log(Si)+log(S2)  (plus  a  constant  to  normalize  to  zero  at  20 
seconds).  As  the  figure  shows,  this  is  a  good  approximation  to  the  average  excitation  function 
across  the  frequency  band,  however  there  is  substantial  regional  variation  in  the  function  that  is 
accounted  for  in  the  path  corrected  magnitudes. 

Figure  2.1  (right)  shows  that  attenuation  calculated  from  earth  models  is  somewhat  higher  than 
the  Rezapour  and  Pearce  attenuation,  and  both  are  higher  than  the  Russell  attenuation,  which  is 
based  on  the  earlier  model  of  von  Seggem  (1977).  The  model-based  attenuation  corresponds  to  a 
Rayleigh  wave  Q  of  about  400,  while  the  Rezapour/Pearce  and  von  Seggem/Russell  attenuation 
correspond  to  Rayleigh  wave  Q  of  about  550  and  800,  respectively.  The  Q  value  of  400  is  more 
consistent  with  empirical  Rayleigh  wave  Q  studies  than  the  higher  values  of  the  other 
magnitudes,  and  later  in  this  report  we  find  that  Q  values  for  much  of  Eurasia  are  still  lower.  The 
higher  Q  values  in  the  Rezapour  and  Pearce  and  von  Seggem  studies  may  be  because  those 
magnitudes  were  based  on  Rayleigh  wave  amplitudes  covering  a  large  distance  range,  and 
Rayleigh  waves  along  lower  Q  paths  have  attenuated  away  at  the  larger  distances,  biasing  the 
attenuation  estimates  to  higher  Q  values. 
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Figure  2.1  Left:  Mean  log(Sl^S2)  and  ±1  standard  deviation  (bine).  Red  line  shows  the  Rnssell  approximation 
for  the  source  function.  Right:  Mean  attenuation  and  ±1  standard  deviation  (bine)  derived  from 
earth  models,  Rnssell  approximation  for  attennation  (red),  and  Rezaponr  and  Pearce  global  estimate 
at  20  seconds  (black).  In  both,  green  marks  show  values  for  individual  Eurasian  structures. 


Figure  2.2  compares  the  magnitude  distance  corrections,  which  are  the  equations  listed  in  Table 

2.1  calculated  with  A=l,  Gmin=0.6,  f  =  and  A  _  fc'^^  ^  20  and  10 

rVA  G 


seconds,  respectively.  Differences  between  the  distance  corrections  are  generally  small.  The 
main  differences  are  the  larger  correction  at  close  distances  for  the  Rezapour  and  Pearce 
magnitude,  and  the  larger  correction  for  the  path  corrected  magnitude  with  model-based 
attenuation  at  large  distances.  Ms(bp)  will,  of  course,  vary  for  each  source  and  receiver  location 
corresponding  to  the  particular  earth  structure  and  path  attenuation.  The  magnitude  correction  at 
close  distances  is  also  larger  for  Ms(b)  than  for  Ms(bp)  because  the  difference  in  attenuation  causes 
a  small  difference  in  the  normalization  constant  which  is  calculated  at  50  degrees. 


Surfece  Wave  Magnilitde  Correction  {20  secorids) 


- Mstbp){MotJel-baseci  gamma) 

- Ms(bp){R(P  gamma) _ 

'^■^0  20  40  eo  eo  lOO 

Distance  (degrees) 


Distance  (degrees) 


Figure  2.2.  20  second  magnitude  correction  vs.  distance  for  Rezapour/Pearce  (dashed  blue),  Ms(b)  (dot-dash 
black),  Ms(bp)  with  model-based  gamma  (solid  red),  and  Ms(bp)  with  Rezapour  and  Pearce  gamma 
(solid  maroon)  (left).  10  second  magnitude  correction  vs.  distance  Ms(b)  (dot-dash  black),  and  Ms(bp) 
with  model-based  gamma  (solid  red)  (right). 
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2.2  Analysis  of  surface  waves  from  the  North  Korean  nuclear  test 

To  examine  the  differences  between  magnitudes  in  more  detail  and  illustrate  some  potential 
problems,  we  apply  the  magnitude  methods  to  surface  waves  from  the  North  Korean  nuclear  test 
of  October  9,  2006.  This  is  a  good  test  case  because  the  surface  waves  are  small  -  above  noise 
level  at  only  7  of  the  closest  stations  (Figure  2.3),  difficult  to  see  at  all  in  the  unfiltered  records, 
but  visible  at  these  stations  when  low  pass  filtered  (Figure  2.4). 


Figure  2.3.  Location  of  the  North  Korean  nuclear  test  and  recording  stations. 
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Figure  2.4.  Data  from  the  North  Korean  explosion  filtered  from  0.01-0.1  Hz.  Surface  waves  are  clearly  visible 
at  all  stations.  HIA  has  a  glitch  or  interfering  arrival  after  the  explosion  arrival.  The  explosion 
arrival  is  visible  just  after  the  BJT  arrival. 


Figure  2.5  shows  the  predicted  spectra  at  each  of  the  7  stations  based  on  the  model  based  source 
and  path  corrections  -  these  are  the  negative  of  the  sum  of  corrections  in  row  5  of  Table  2.1,  and 
are  equivalent  to  predicted,  normalized  explosion-generated  surface  wave  spectra  at  each 
location.  Differences  between  the  model-based  and  Russell  sets  of  corrections  (Figure  2.5,  right) 
range  from  -0.15  to  0.05  magnitude  units. 

Path  Corrected  Magnitude  Can'ectiqn  Path  Corrected  Magnitude  Correcticn  Difference 


Figure  2.5.  Path  corrections  for  the  7  stations  that  recorded  surface  waves  from  the  North  Korean  nuclear 
test  using  model  based  corrections  (left)  and  differences  between  the  model-based  and  Russell  path 
corrections  for  the  7  stations  recording  the  North  Korean  event. 

Figure  2.6  (left)  shows  the  calculated  Butterworth  filtered  and  path  corrected  spectral  magnitudes 
for  station  BJT  for  6  values  of  filter  width,  as  specified  by  Gmin  ranging  from  0. 1  to  0.6.  The 
higher  values  correspond  to  wider  filter  widths  which  has  the  effect  of  smoothing  the  spectrum 
and  giving  more  consistent  values  between  frequencies.  A  disadvantage  of  the  larger  values. 
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however,  is  that  the  frequency  band  extends  farther  outside  of  the  band  of  interest,  possibly 
allowing  contamination  by  noise  or  other  phases.  Figure  2.6  (right)  shows  the  Butterworth 
fdtered  magnitude  for  Gmin=0.6,  the  path  corrected  magnitude,  the  path  corrected  spectral 
magnitude,  and  the  best  value  of  the  path  corrected  spectral  magnitude  calculated  using  a  robust 
mean  (Stevens  et  al.,  2005).  Figure  2.7  shows  a  comparison  of  Butterworth  fdtered  magnitudes 
with  and  without  path  corrections.  If  the  procedure  were  working  perfectly  and  the  surface  wave 
spectra  were  just  like  synthetics,  then  all  of  the  curves  in  Figures  2.6  and  2.7  would  be  flat  lines. 
The  path  corrected  spectra  are  slightly  flatter  than  the  Butterworth  fdtered  spectra,  but  it  is  clear 
that  unmodeled  variations  are  significantly  larger  than  the  differences  between  the  individual 
magnitude  curves. 


BJT  mspc  filter  widlb  gmin  BJT  Spectral  Magnitude 


Figure  2.6.  Calculation  of  path-corrected  Butterworth  filtered  magnitude  for  6  filter  widths  (left)  at  station 
BJT.  Butterworth  filtered,  path  corrected,  and  spectral  magnitudes  at  BJT  (right). 


msbp.  G=Q.6  mSpc,  G=0.6 


Figure  2.7.  Butterworth  filtered  (left)  and  path  corrected  (right)  magnitudes  for  the  North  Korean  nuclear 
test  with  varying  filter  widths. 

The  spectra  divide  into  two  distinct  groups  -  lower  values  at  INCN,  KSRS  and  MDJ,  and  higher 
values  at  BJT,  ENH  and  TLY.  HIA  is  more  complicated,  but  generally  higher.  Furthermore,  the 
three  lower  amplitude  stations,  which  are  also  the  closest  stations,  have  much  flatter  corrected 
spectra  than  the  three  higher  amplitude  stations,  suggesting  a  possible  frequency  dependent 
amplification  of  these  surface  waves.  Also,  the  lower  amplitude  stations  are  located  due  north 


11 


and  south  of  the  event,  while  the  other  stations  are  located  to  the  west  (and  northwest/southwest) 
Since  we  have  corrected  for  source  and  receiver  structure,  and  attenuation  differences  could  not 
be  responsible  for  differences  this  large  on  paths  this  short,  there  are  two  remaining  likely  causes 
for  the  amplitude  variations:  tectonic  release  (or  other  azimuth  dependent  source  component), 
and  focusing  due  to  path  structure. 

The  effect  of  path  structure  calculated  using  the  Bom  approximation  (Zhou  et  al.,  2004)  and  the 
earth  models  of  Stevens  et  al.  (2005)  is  shown  in  Figure  2.8.  Although  there  are  clearly  some  big 
differences  between  the  predictions  and  observations,  there  are  also  some  interesting  similarities. 
First,  the  corrections  for  the  three  closest  stations  are  almost  flat  and  separated  by  approximately 
the  same  amount  as  the  observations.  Second,  the  amplitude  correction  for  ENH  is  similar  to 
what  is  observed  for  ENH  and  also  similar  to  the  observations  at  TLY  and  BIT.  That  is,  the  data 
show  a  peak  in  the  spectmm  in  the  middle  periods,  dropping  back  close  to  the  level  of  the  closer 
stations  at  the  longest  and  shortest  periods.  On  the  other  hand  a  very  large  amplification,  which  is 
not  observed,  is  predicted  for  BJT  due  to  a  grazing  path  along  the  north  end  of  the  Yellow  Sea, 
and  a  big  decrease,  also  not  observed,  is  predicted  for  TLY.  There  are  two  likely  explanations  for 
this:  1)  the  one  degree  stmctural  resolution  is  not  sufficient  for  these  short  paths;  and/or  2)  the 
stractural  complexity  exceeds  the  limits  of  the  Bom  approximation.  The  amplification  predicted 
for  BJT  by  the  Bom  approximation,  for  example,  is  very  strongly  dependent  on  exactly  where 
the  station  is  with  respect  to  the  low  velocity  zone.  The  Bom  approximation  and  an  analysis  of 
the  effectiveness  of  its  amplitude  correction  are  discussed  in  detail  in  the  next  section  of  this 
report. 


Predicted  Born  Amplitude  Enhancement 


Figure  2.8.  Amplitude  corrections  due  to  structure  predicted  using  the  Born  approximation. 

Tectonic  release  could  cause  the  observed  ^factor  of  2  offset  between  the  two  station  sets.  If,  for 
example,  tectonic  release  had  the  mechanism  of  a  normal  fault  with  tension  in  the  east-west 
direction,  it  would  amplify  stations  to  the  west  more  than  to  the  north.  However,  to  cause  a  factor 
of  2  difference  would  require  a  secondary  source  large  enough  that  large  Love  waves  would  also 
be  expected.  No  Love  waves  are  apparent  in  any  of  the  data.  In  any  case,  path  stmcture  seems  a 
more  likely  cause  of  the  variability,  and  it  is  possible  that  the  variations  would  be  predictable 
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with  higher  resolution  earth  models,  and/or  a  more  accurate  method  for  accounting  for  structural 
variations. 

One  final  issue  is  the  most  appropriate  measurement  of  Ms  derived  from  spectra.  In  our  previous 
work,  and  in  the  discussion  above  of  the  path  corrected  spectral  magnitude,  we  have  used  a 
robust  mean  value  as  the  best  estimate.  Bonner  et  ah,  however,  have  instead  used  the  maximum 
value  of  the  Butterworth  filtered  Mj.  This  has  some  advantage  in  discrimination,  because 
although  the  corrections  flatten  explosion  spectra,  they  do  not  quite  flatten  earthquake  spectra, 
which  tend  to  be  larger  at  lower  frequencies  (Figure  2.9).  So  in  principle  the  Mg  value  for 
explosions  is  frequency  independent,  while  the  earthquake  Mg  can  choose  a  higher  value, 
improving  the  Mgimb  discriminant.  As  the  explosion  spectra  above  demonstrate,  however,  there 
are  variations  in  the  explosion  spectra  also,  and  choosing  the  maximum  point  can  lead  to  an  Mg 
that  is  too  high,  degrading  discrimination.  We  made  one  modification  to  Bonner  et  al’s 
Ms(VMAX)  procedure  by  implementing  a  simple  outlier  rejection  test  -  we  calculate  the  mean 
and  standard  deviation  of  the  magnitudes  and  reject  outliers  greater  than  two  standard  deviations 
from  the  mean.  We  then  recalculate  the  mean  and  repeat  the  outlier  rejection,  using  the 
remaining  points  to  calculate  either  the  mean  or  maximum  magnitude.  Tables  2.2  and  2.3  show  a 
comparison  of  magnitudes  for  this  event  using  the  different  methods  (defined  in  Table  2.1).  The 
results  show  that  the  mean  value  estimates  are  more  consistent  than  the  peak  value  estimates. 

Table  2.2.  Station  and  network  mean  valnes  for  M,  calculated  nsing  the  mean  values  over  the  frequency  band. 


Station 

Distance  (km) 

Ms(b)  Mean 

Ms(bp)  Mean 

Ms(sp)  Mean 

MDJ 

369 

2.65 

2.67 

2.69 

KSRS 

440 

2.69 

2.70 

2.76 

INCN 

476 

2.68 

2.69 

2.81 

BJT 

1103 

3.03 

3.02 

3.07 

HIA 

1148 

2.93 

2.93 

3.00 

ENH 

2147 

3.01 

3.06 

3.20 

TLY 

2252 

3.00 

3.09 

3.16 

Mean  (SD) 

2.86  (0.18) 

2.88  (0.19) 

2.96  (0.20) 

Table  2.3.  Station  and  network  mean  valnes  for  M,  calculated  using  peak  valnes  within  the  frequency  band. 


Station 

Distance  (km) 

Ms(b)  Peak 

Ms(bp)  Peak 

MDJ 

369 

2.71 

2.73 

KSRS 

440 

2.75 

2.74 

INCN 

476 

2.80 

2.78 

BJT 

1103 

3.18 

3.13 

HIA 

1148 

3.07 

3.08 

ENH 

2147 

3.21 

3.23 

TLY 

2252 

3.25 

3.32 

Mean  (SD) 

3.00  (0.23) 

3.00  (0.25) 
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Scalar  Moment  Estimate 


Frequency  (FIz) 

Figure  2.9.  Path  corrected  spectra  for  an  explosion  and  for  earthquakes  calculated  for  several  depths.  The 
path  corrected  explosion  spectrum  is  flat  over  the  entire  frequency  band  (for  perfect  data  and  path 
correction),  while  the  path  corrected  earthquake  spectrum  is  flattened,  but  has  some  variation  due  to 
source  mechanism  and  source  depth  and  generally  decreases  with  increasing  frequency. 

2.3  Yield  estimation,  earthquake/explosion  discrimination  and  the  North  Korean  test 

The  Ms  measurements  for  this  event  are  larger  than  expected  for  an  explosion  with  mb  4.1.  Based 
on  typical  mbiyield  curves  for  hard  rock,  an  event  of  that  size  would  have  a  yield  of  less  than  one 
kiloton.  However,  Stevens  and  Murphy  (2001)  show  that  Mjiyield  data  are  fit  well  with  a 
relation  Ms=log(Y)+2.1,  which  for  Ms  2.9  gives  an  estimated  yield  of  about  6  kilotons.  On  the 
other  hand,  the  Mjimb  curve  for  granite  from  Stevens  and  Day  (1985)  comes  quite  close  to  the 
Ms  and  mb  for  the  North  Korean  event.  The  yield  estimate  for  this  event  and  corresponding 
magnitudes  are  also  affected  by  the  relatively  deep  source  depth.  Bennett  et  al  (2006)  estimate 
the  source  depth  between  200  and  800  meters  based  on  the  turmel  location  and  topography.  They 
estimate  a  yield  of  about  one  kiloton  using  a  comparison  of  P-wave  spectra  with  a  source  depth 
of  800  meters.  In  this  section  we  examine  the  differences  between  these  estimates. 

Figure  2.10  shows  Mjimb  data  from  historical  explosions  together  with  theoretical  curves  for  Ms 
and  mb  from  Stevens  and  Day  (1985)  for  Mueller-Murphy  granite  (Murphy,  1977).  Curves  are 
shown  for  normal  containment  depth  of  122Y  meters  and  for  the  estimated  800  meter  depth  of 
this  event.  Although  the  Korean  data  point  appears  to  agree  quite  well  with  the  estimate,  this  is 
somewhat  misleading  because  the  Stevens  and  Day  magnitudes  were  normalized  to  an  NTS 
magnitude  yield  curve,  and  so  the  mb  estimates  are  biased  low  compared  to  Eurasian  mb  values. 
Figure  2.10  also  shows  an  Mjimb  curve  shifted  by  0.4  in  mb,  which  is  more  appropriate  for  the 
North  Korean  location.  So  the  Ms  for  this  event  is  still  high  relative  to  predictions,  but  consistent 
with  the  upper  part  of  the  distribution  of  other  explosions. 
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M  vs, 

s  b 


Figure  2.10.  Msinib  curves  for  Mueller  Murphy  granite  using  formulas  from  Stevens  and  Day  (1985)  and  with 
0.4  mb  bias  correction.  Data  from  Stevens  and  Murphy  (2001),  plus  the  North  Korea  data  point. 

The  nominal  mbiyield  relation  for  Eurasian  hardrock  is  mb=0.75  logioY  +  4.45  (Murphy,  1996), 
which  would  imply  a  yield  of  less  than  /a  kiloton  for  this  event.  Figure  2.1 1  shows  a  comparison 
of  the  nominal  mbiyield  relation  and  the  granite  Mueller-Murphy  mbiyield  curve  using  standard 
depth  and  the  800  meter  depth  with  the  bias  adjusted  Stevens  and  Day  (1985)  relations.  One 
kiloton  is  in  good  agreement  with  these  curves  for  overburied  depths. 
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vs.  Yield  at  BOO  meters 

b 


Figure  2.11.  lUbiyield  curves  for  Mueller  Murphy  granite  at  normal  containment  depth  and  at  800  m  depth 
using  formulas  from  Stevens  and  Day  (1985),  and  nominal  Eurasian  mbiyield  curve  (Murphy,  1996). 
Data  from  Stevens  and  Murphy  (2001),  plus  North  Korea  mi,  and  yield  estimate. 


M  vs.  Yield 
s 


Figure  2.12.  Msiyield  curves  for  Mueller  Murphy  granite  using  formulas  from  Stevens  and  Day  (1985)  for 
normal  containment  depth  and  800  meter  depth,  and  Msiyield  curve  and  data  from  Stevens  and 
Murphy  (2001),  plus  North  Korea  Ms  and  yield  estimate. 

Figure  2.12  shows  three  Ms: Yield  curves:  the  best  fit  Mj: Yield  curve  with  unit  slope  from 
Stevens  and  Murphy  (2001),  and  Ms:Yield  using  Mueller-Murphy  granite  and  the  Stevens  and 
Day  (1985)  relations  with  normal  scaled  depth  and  a  depth  of  800  meters.  Although  the  Ms 
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appears  to  agree  well  with  the  normal  containment  depth  curve,  as  discussed  above  it  is  very 
likely  that  the  explosion  was  overburied.  The  main  difference  between  the  three  curves  in  Figure 
2.12  is  the  slope  of  the  lines.  The  Stevens  and  Murphy  curve  assumed  a  unit  slope,  which  is 
expected  for  varying  yield  at  fixed  depth.  The  Mueller-Murphy  model  has  an  empirically 
determined  slope  of  0.87  at  fixed  depth,  and  a  slope  of  0.76  when  normal  containment  depth  is 
assumed.  The  smaller  slope  occurs  because  lower  overburden  pressure  at  shallower  depths 
increases  the  source  size. 

While  the  North  Korean  data  point  is  more  consistent  with  the  Mueller-Murphy  curves,  the  Ms  is 
still  high  by  about  0.4  magnitude  units  relative  to  the  Mueller-Murphy  curve  for  an  overburied 
explosion,  and  it  isn’t  clear  whether  the  North  Korean  data  point  indicates  a  slope  less  than  one 
in  general,  or  whether  this  particular  data  point  is  at  the  high  end  of  the  MjiYield  distribution. 
There  are  two  factors  that  could  contribute  to  the  apparently  high  Ms  relative  to  the  historical 
data  set.  The  historical  data  is  dominated  by  events  at  NTS  and  Balapan  and  so  1)  the  lower 
material  velocities  at  NTS  increase  the  Msimb  separation  as  discussed  by  Stevens  and  Day,  1985; 
and  2)  many  of  the  events  at  Balapan  exhibit  compressive  tectonic  release,  reducing  Ms  for  those 
events  (Day  and  Stevens,  1986).  The  combination  of  these  two  effects  causes  Ms  for  a  small 
event  in  hard  rock  with  no  tectonic  release  to  be  higher  than  predicted  from  a  curve  with  slope  1 
derived  from  a  data  set  with  a  mixture  of  source  media  velocities  and  a  substantial  number  of 
events  with  Ms  reduced  by  tectonic  release. 


17 


3.0  Amplitude  Corrections  Using  the  Born  Approximation 

Surface  wave  amplitudes  are  affected  by  both  attenuation  and  earth  structure.  The  effect  on 
surface  wave  amplitudes  of  propagation  normal  to  variations  in  earth  structure  is  predicted  fairly 
well  by  conservation  of  energy.  Propagation  along  paths  at  grazing  incidence  to  large  structure 
variations,  however,  are  much  more  difficult  to  predict.  Our  main  interest  in  this  project  is  on 
understanding  amplitude  variations  in  8-15  second  surface  waves.  In  this  frequency  band,  surface 
waves  may  be  affected  as  strongly  or  more  strongly  by  earth  structure  than  by  intrinsic 
attenuation,  particularly  along  shorter  paths.  Our  goal  is  therefore  to  be  able  to  model  and 
correct  for  both  of  these  effects.  Our  approach  is  illustrated  in  Figure  3.1. 


Figure  3.1.  Overview  of  surface  wave  dispersion  and  attenuation  analysis.  We  derive  earth  structure  from  a 
large  number  of  dispersion  curves,  then  use  the  Born  approximation  to  derive  amplitude  corrections 
using  this  structure  prior  to  Q  inversion.  The  Q  structures  are  then  used  to  predict  attenuation  along 
any  path. 

In  an  earlier  project  (Stevens  et  ah,  2005)  we  developed  global,  regionalized  dispersion  models 
that  allow  the  phase  and  group  velocity  to  be  calculated  between  any  two  points  on  the  earth.  We 
did  this  by  accumulating  a  large  data  set  consisting  of  more  than  1  million  dispersion 
measurements  derived  by  a  number  of  researchers,  and  then  inverting  this  data  set  to  determine 
earth  structure,  which  in  turn  was  used  to  generate  dispersion  maps  at  all  frequencies.  In  that 
project,  we  modeled  surface  waves  in  a  heterogeneous  earth  using  the  following  approximations: 
1)  surface  waves  propagate  along  great  circle  paths,  2)  surface  wave  phase  and  group  velocities 
and  anelastic  attenuation  can  be  modeled  using  a  path  integral  between  source  and  receiver,  and 
3)  energy  is  conserved  with  no  mode  conversion  across  material  boundaries.  This  approximation 
is  quite  good  for  large  parts  of  the  world,  particularly  at  lower  frequencies,  but  the  unmodeled 
variations  become  important  in  regions  of  structural  complexity. 

3.1  Surface  wave  amplitude  predictability 

An  important  goal  of  this  project  is  to  be  able  to  predict  surface  wave  amplitudes  in  both  simple 
and  complex  structures,  to  determine  under  what  conditions  the  more  complicated  calculations 
required  for  laterally  heterogeneous  structure  are  required,  and  under  what  conditions  the 
approximations  generally  used  for  calculating  surface  waves  in  complex  structures  become 
inadequate.  In  the  following,  we  discuss  calculations  of  surface  waves  in  simple  and  complex 
structures. 
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3.2  Surface  wave  propagation  in  simple  structures 

We  define  “simple  structures”  to  mean  those  structures  in  which  the  surface  wave  propagation  is 
normal  to  all  changes  in  structure,  and  lateral  changes  in  structure  are  negligible.  In  that  case  we 
can  predict  surface  wave  amplitude  and  phase  using  an  approximation  originally  due  to  McGarr 
(1969)  that  uses  propagation  of  surface  waves  along  great  circle  paths  with  conservation  of 
energy  across  material  interfaces  and  no  mode  conversion.  With  these  approximations,  surface 
wave  propagation  in  a  heterogeneous,  anelastic  structure  takes  the  following  form,  separating 
source,  path  and  receiver  (notation  follows  Harkrider  et  al,  1994): 


u^{co,r,(p) 


Yp^ 


(3.1) 


where  co  is  angular  frequency,  r  is  source  to  receiver  distance,  h  is  source  depth,  ae  is  the  radius 
of  the  earth,  (p  is  azimuth,  is  the  Rayleigh  wave  amplitude  function,  c  is  phase  velocity,  y  is 
the  attenuation  coefficient,  and  the  subscripts  1,  2,  and  p  refer  to  parameters  derived  from  the 
source  region  structure,  parameters  derived  from  the  receiver  region  structure,  and  parameters 
which  are  defined  by  path  averages,  respectively.  All  source  properties  are  contained  in  the 
function  Fs.  For  an  isotropic  explosion  source,  the  Rayleigh  wave  spectrum  can  be  written: 

w,  {co,r)  =  M, -  '  ,  , - - -  (3  -2) 

yja^sm(r/  aj 


where  (po  is  the  initial  phase  equal  to  -37i/4,  depends  on  the  source  region  elastic  structure  and 
the  explosion  source  depth,  depends  on  the  receiver  region  elastic  structure.  where 

Mo  is  the  explosion  isotropic  moment.  is  defined  this  way  so  that  the  function  does  not 

depend  explicitly  on  the  material  properties  at  the  source  depth.  (More  details  are  given  in 
Stevens  and  McLaughlin  (2001)  and  Stevens  and  Murphy  (2001)). 

3.3  Surface  wave  propagation  in  complex  structures 

In  more  complex  structures,  equations  3.1  and  3.2  may  be  inadequate  to  describe  surface  waves. 
In  our  experience,  these  approximations  work  very  well  for  modeling  surface  wave  dispersion 
and  fairly  well  for  modeling  surface  wave  amplitude  over  most  of  the  world.  In  this  project, 
however,  we  are  performing  a  more  complete  analysis  including  an  approximation  for  the  effects 
of  scattering  and  diffraction.  This  is  important  for  two  reasons: 

1)  Some  of  the  remaining  residual  in  the  global  dispersion  models  is  due  to  scattering  and 
diffraction,  and  incorporation  of  these  effects  into  our  analysis  may  allow  us  to  correct 
for  them;  and 

2)  To  perform  inversion  of  attenuation  data  for  Q  structure  as  described  in  the  following 
section,  we  want  to  correct  the  amplitude  for  the  effects  of  heterogeneous  structure.  The 
effect  of  heterogeneous  structure  on  amplitude  is  stronger  than  on  dispersion. 

Modeling  of  scattering  and  diffraction  is  an  active  area  of  current  research.  Most  of  the  research 
relevant  to  this  project  use  variants  of  the  single-scattering  Bom  approximation  to  model  the 
scattered  wave  field  (Snieder,  1986).  Zhou  et  al.  (2004)  summarize  this  work  and  derive 
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sensitivity  kernels  for  phase,  amplitude  and  arrival  angle.  The  Bom  approximation  models  the 
observed  surface  wave  at  a  receiver  as  a  sum  of  a  direct  wave  plus  waves  scattered  from  material 
inhomogeneities  throughout  the  region. 

Ritzwoller  et  al.  (2002)  used  a  simplified  version  of  the  Bom  approximation  to  include 
diffraction  in  surface  wave  tomography.  They  modeled  the  sensitivity  kernel  with  a  boxcar 
function  the  width  of  the  Fresnel  zone  normal  to  the  source  to  receiver  path,  then  used  an  area 
integral  over  this  region  in  place  of  the  ray  theory  path  integral  for  performing  tomographic 
inversion.  Yoshizawa  and  Kennett  (2004)  and  Kennett  and  Yoshizawa  (2002)  use  a  similar 
technique  with  a  narrower  kernel  that  they  believe  to  be  more  representative  of  realistic  surface 
waves.  Spetzler  et  al.  (2001,  2002)  discuss  the  implications  of  scattering  and  diffraction  for 
surface  wave  tomography.  Friedrich  et  al.  (1993),  Friedrich  (1999)  and  Maupin  (2001)  extend 
the  Bom  approximation  to  incorporate  multiple  scattering. 


We  implemented  the  algorithms  of  Zhou  et  al.  (2004)  for  calculating  finite  frequency  sensitivity 
kernels  for  dispersion  and  amplitude  variations.  Using  the  forward  scattering,  forward 
propagating  approximation,  the  phase  and  amplitude  corrections  are: 


||x'^(Jc/c)(iQ,  where 


Ik^'^  sin  [/t  ( A '+  A A)  +  ^4] 
^8;r  sin  A '  sin  A "/ sin  A 


(3.3) 


^In  ^  {5clc)dQ. ,  where 


cos \k  (A  ’+  A A)  +  V4] 
^8 sin  A '  sin  A "/ sin  A 


(3.4) 


where  distance  is  in  radians,  k  is  wavenumber,  and  A',  A"  and  A  refer  to  the  source  to  scatterer, 
scatterer  to  receiver,  and  source  to  receiver  distances,  respectively.  The  integrals  ran  over  the 
entire  earth’s  surface,  although  in  practice  (and  in  this  report)  are  limited  to  the  first  Fresnel 
zone,  which  is  defined  by  A:(A'  +  A"  -  A)  <  37r/4.  Phase  and  amplitude  kemals  at  12  seconds  are 
shown  in  Figure  3.2.  Dahlen  and  Zhou  (2006)  extend  these  equations  to  derive  group  delay  and 
intrinsic  attenuation  kernels. 


Figure  3.2. 12  second  phase  velocity  (left)  and  amplitude  (right)  kernels  across  central  Asia. 
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3.4  Application  of  corrections  to  surface  wave  amplitudes 

We  tested  the  algorithms  described  above  using  the  one  degree  dispersion  maps  of  Stevens  et  al. 
(2005).  While  the  results  are  reasonable  for  prediction  of  dispersion  variations,  the  predicted 
amplitude  corrections  seem  unreasonably  large,  particularly  on  long  paths.  Consequently,  we 
have  performed  some  additional  calculations  and  analysis  to  investigate  how  model  roughness 
affects  amplitudes.  Figure  3.3  shows  the  10  second  phase  velocity  model  for  Eurasia;  Figure  3.4 
shows  the  predicted  amplitude  variation  for  paths  through  this  region  from  the  Lop  Nor  test  site 
using  the  model  shown  in  Figure  3.3  directly,  and  using  a  “smoothed”  version  of  the  amplitude 
variation  in  which  the  phase  slowness  was  modeled  with  a  bilinear  function  instead  of  discrete 
blocks.  For  both  of  the  amplitude  figures,  the  anomalies  have  been  truncated  where  they  exceed 
logio(amplitude)  =  0.6.  The  amplitude  variations  become  quite  large  on  longer  paths,  and  it 
does  not  appear  that  the  Bom  approximation  is  giving  reasonable  answers  on  these  paths. 


10  s  Phase  Velocity 


Figure  3.3.  Eurasian  phase  velocity  model  at  10  seconds  from  Stevens  et  al.  (2005). 
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Figure  3.4.  Left  -  predicted  amplitude  variations  at  10  seconds  through  the  phase  velocity  model  of  Figure  2 
on  paths  out  of  Lop  Nor.  Right  -  same,  but  the  velocity  model  has  been  smoothed  by  modeling  it  as  a 
bilinear  rather  than  piecewise  discontinuous  function. 
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In  order  to  validate/correct  the  results,  we  need  to  determine  whether  the  Bom  approximation  is 
giving  correct  results  for  stractural  variations  of  the  magnitude  present  in  the  Stevens  et  al. 
(2005)  models,  and  at  the  frequencies  of  interest.  To  do  this,  we  performed  a  test  case  of  a 
stmcture  for  the  Tarim  Basin  embedded  in  a  uniform  stmcture  typical  of  shield  regions  of 
Eurasia,  such  as  those  that  surround  the  Tarim  Basin  (Figure  3.5).  We  performed  a  Born 
approximation  calculation  and  a  3D  finite  difference  calculation  for  a  source  located  just  east  of 
the  Tarim  Basin  and  examined  the  wavefield  for  several  hundred  km  west  of  the  Tarim  Basin. 
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Figure  3.5.  Comparison  of  Born  approximation  (left)  with  finite  difference  calculation  (right)  of  amplitude 
perturbations  at  20  seconds  (top)  and  10  seconds  (bottom).  The  rectangular  inclusion  is  modeled 
after  the  Tarim  Basin  structure,  and  the  external  structure  after  a  Eurasian  shield  earth  structure. 
The  source  is  on  the  horizontal  axis  at  the  right  edge  of  the  plot.  There  is  general  agreement  in  the 
features  of  the  two  calculations.  The  amplitude  is  increased  in  a  band  above  and  to  the  left  of  the 
inclusion  in  both  cases,  and  decreased  above  that.  However,  there  are  some  interference  effects  in  the 
finite  difference  calculation  that  are  not  reproduced  in  the  Born  calculation,  and  the  patterns  of  high 
and  low  amplitudes  to  not  match  exactly. 
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The  results  showed  that  although  the  Bom  approximation  is  generally  consistent  with  the  finite 
difference  results,  there  are  noticeable  interference  effects  leading  to  high  and  low  amplitudes  in 
the  finite  difference  calculation  that  are  not  present  in  the  Bom  approximation.  The  differences 
are  significantly  larger  at  10  seconds  than  at  20  seconds. 

One  reason  for  this  increase  in  complexity  is  illustrated  in  Figure  3.6.  Propagation  of  the 
cylindrical  wave  leaving  the  source  through  the  Tarim  Basin  model  leads  to  a  strong  diffracted 
wave  generated  by  the  wavefront  passing  along  the  top  of  the  Basin.  This  secondary  wave 
interferes  with  the  direct  wave  and  complicates  analysis,  particularly  in  the  interpretation  of 
spectra.  Since  the  first  order  Bom  approximation  only  models  the  direct  wave,  it  cannot 
reproduce  this  strongly  diffracted  secondary  wave,  although  it  may  do  an  adequate  job  of 
predicting  the  primary  arrival.  Also  shown  in  Figure  3.6  are  two  observed  waveforms  that 
traveled  through  the  Tarim  Basin.  There  are  two  distinct  surface  wave  arrivals  similar  to  the 
figure  on  the  left.  Although  we  have  not  done  sufficient  analysis  to  say  that  the  split  in  these 
seismograms  was  due  to  the  effect  illustrated  in  the  left  figure,  it  does  suggest  that  strong 
diffraction  may  be  responsible. 


0  3»  4«  e»  tooo  iaw 


Figure  3.6.  Vertical  component  velocity  after  propagation  across  the  low  velocity  basin  (left).  There  is  a 
strong  diffracted  wave  that  interferes  with  the  direct  wave.  The  right  figure  shows  two  observed 
waveforms  that  passed  through  the  Tarim  Basin  and  have  two  distinct  surface  wave  arrivals. 
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3.4  Analysis  of  surface  waves  at  KNET  propagating  through  and  near  the  Tarim  Basin 

The  Tarim  Basin  is  a  particularly  troublesome  region  for  surface  wave  propagation  because  of 
the  presence  of  deep  sediments  with  variable  thickness  (Figure  3.7).  We  examined  surface  waves 
in  this  region  specifically  to  see  if  the  kind  of  amplitude  variations  predicted  by  the  Bom 
approximation  are  observable  here.  Figure  3.7  shows  the  paths  between  events  southeast  of  the 
Tarim  Basin  that  were  recorded  at  KNET.  Figure  3.8  shows  the  seismograms  recorded  along 
these  paths  over  three  frequency  bands,  and  Figure  3.9  shows  the  seismograms  from  the  same 
events  recorded  along  a  northerly  path  that  is  much  less  complex.  The  seismograms  do  in  fact 
show  what  we  would  expect:  the  paths  through  and  near  the  Tarim  Basin  show  much  more 
complexity  than  the  seismograms  along  the  simpler  paths.  There  is  more  variability  at  high 
frequency  (10-20  seconds)  and  a  strong  amplification  along  the  southernmost  grazing  path, 
similar  to  what  we  would  expect  from  a  seismogram  in  the  red  bands  of  amplification  in  Figure 
3.5.  So  the  seismograms  are  qualitatively  similar  to  what  is  predicted  by  the  Bom  approximation 
(or  finite  difference  calculations)  for  a  complex  region.  In  the  next  section  we  examine  the 
predictability  of  these  amplitude  variations. 


Figure  3.7.  Events  within  green  circle  (lower  right)  propagate  at  various  distances  in  from  the  southern  edge 
of  the  Tarim  Basin  to  a  KNET  station  (red  lines).  The  northernmost  event  proved  to  have  a  complex 
source.  All  other  events  are  shown,  ordered  by  azimuth  (northernmost  path  topmost),  in  the  next  two 
images.  Scale  indicates  sediment  thickness  in  km. 


24 


40  -  50  Seconds  20  -  30  Seconds  10-20  Seconds 


0  2  4  6  0  2  4  6  0  2  4 

X  10+2  X  10+2  X  10+2 


Figure  3.8  KNET  seismograms  of  the  events  with  paths  through  the  Tarim  Basin.  All  records  are  normalized 
by  the  40-50  second  period  surface  wave  amplitude.  Amplitudes  vary  at  higher  frequency,  as 
predicted  by  the  simulations.  The  amplitude  of  the  10-20  second  period  surface  wave  from  the  event 
with  the  most  glancing  path  becomes  very  large  (bottom  trace). 
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Figure  3.9.  Seismograms  from  the  same  events  as  the  previous  figure,  in  the  same  order  and  normalized  by 
the  40-50  second  surface  wave  amplitudes,  but  recorded  at  KMI  (-130”  azimuth,  away  from  the 
Tarim  Basin).  The  amplitudes  are  more  consistent  and  waveforms  much  less  complex,  indicating  that 
the  amplitude  variation  and  complexity  observed  at  KNET  can  be  attributed  to  the  effects  of  the 
propagation  near  the  basin  boundary. 


25 


3.5  Analysis  of  surface  waves  at  KNET  stations  due  to  Lop  Nor  explosions 

To  assess  the  predictability  of  amplitude  variations  using  the  Bom  approximation,  we  calculated 
surface  waves  from  several  Lop  Nor  nuclear  tests  recorded  at  KNET  with  and  without  the  Bom 
approximation,  and  then  also  performed  a  large  3D  finite  difference  calculation  of  the  path 
between  Lop  Nor  and  KNET  for  comparison.  The  results  for  all  of  the  Lop  Nor  tests  are  similar 
and  we  show  one  example  in  detail. 

Figure  3.10  shows  seismograms  from  the  mb  5.94  June  8,  1996  explosion  filtered  from  0.01  to 
0.1  Hz  together  with  the  paths  taken  between  Lop  Nor  and  KNET.  Figure  3.11  shows 
uncorrected  and  path  corrected  spectra  (see  section  2)  together  with  the  predicted  Bom 
corrections  for  each  path  calculated  using  the  earth  models  of  Stevens  et  al.  (2005).  The  path 
corrections  do  a  fairly  good  job  of  flattening  the  spectra,  and  the  Bom  corrections  predict 
substantial  variations  in  amplitude  at  higher  frequencies  that  are  not  apparent  in  the  data  (note 
that  the  scales  are  different,  however  -  the  maximum  Bom  correction  is  about  a  factor  of  2, 
although  since  the  corrections  are  additive,  the  negative  corrections  would  be  a  more  substantial 
change).  Figure  3.12  shows  synthetic  seismograms  calculated  for  the  same  paths  with  and 
without  Bom  corrections.  Although  the  Bom  corrections  cause  amplitude  variations  comparable 
to  those  observed,  they  do  not  cause  them  at  the  same  stations.  USP  for  example  is  predicted  to 
be  reduced  substantially  while  the  observed  amplitude  is  relatively  large.  UCH  is  predicted  to  be 
large  and  observed  to  be  small.  Furthermore,  none  of  the  synthetics  match  the  complexity  of  the 
observed  waveforms.  The  predictive  capability  of  the  Bom  approximation  does  not  appear  to  be 
very  good  for  these  complex  paths. 


Figure  3.10.  Paths  (top)  and  data  (bottom)  for  the  Lop  Nor  explosion  of  June  8, 1996.  Data  is  filtered  from 
0.01  to  0.1  Hz. 
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Synthetics  «  Spectra 


igure  3.11.  Raw  (upper  left)  and  Path  Corrected  (upper  right)  Spectra  for  the  Lop  Nor  explosion  of  June  8, 


1996,  and  Born  corrections  (bottom)  calculated  for  the  same  paths.  Perfectly  corrected  spectra 


should  be  flat. 
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Figure  3.12.  Synthetic  surface  wave  seismograms  without  (left)  and  with  (right)  Born  corrections. 
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3.6.  3D  Calculations  of  Lop  Nor  to  KNET  and  comparison  with  the  Born  approximation 

We  use  the  3D  elastic  finite  difference  code  TRES3D  to  calculate  wave  propagation  in  an  earth 
model  corresponding  to  the  region  between  Lop  Nor  and  KNET.  We  measured  surface  wave 
spectral  amplitudes  from  the  calculations  using  the  same  techniques  used  to  measure  observed 
surface  waves:  narrow-band  filtering  to  construct  a  phase-matched  filter  followed  by  phase- 
matched  filtering  to  isolate  the  surface  wave  (Stevens  and  McLaughlin,  2001).  The  numerical 
grid  size  is  678  (North)  X  858  (East)  X  202  (Depth)  covering  an  area  of  approximately  (36-48N) 
and  (71.5-92.5E).  The  grid  spacing  in  each  dimension  is  2km  and  time  step  is  0.12s. 

The  heterogeneous  grid  model  is  generated  from  the  1”  x  1“  global  model  of  Stevens  et  al.  (2005) 
and  the  homogeneous  model  is  obtained  at  the  LopNor  (41.55N,88.70E)  source  site.  No 
attenuation  is  included  in  the  numerical  calculations.  The  homogenous  model  is  shown  in  Figure 
3.13  and  a  slice  through  the  heterogeneous  model  is  shown  in  Figure  3.14. 
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Figure  3.13.  Homogeneous  layered  model.  Numbers  show  Vp,  Vs  and  density.  Trailing  zeros  indicate  inflnite 
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Figure  3.14.  Heterogeneous  model  vertical  profile  of  Vp  (left)  and  Vs  (right)  at  about  42N. 
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The  calculated  amplitude  variations  for  an  explosion  at  Lop  Nor  are  shown  in  Figure  3.15  (20s) 
and  Figure  3.16  (10s).  The  perturbation  for  the  numerical  calculations  is  ln(A/Ao),  A  is  the 
heterogeneous  solution  and  Ao  is  the  homogeneous  solution.  Both  the  numerical  predictions  and 
the  theoretical  predictions  show  complex  variation  patterns  in  comparison  with  those  in  a  simple 
Tarim  Basin  model. 
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Figure  3.15.  Amplitude  perturbations  ln(A/Ao)  at  20  seconds  for  finite  difference  calculation  (left)  and  Born 
approximation  (right).  The  star  is  at  the  source  location. 
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Figure  3.16.  Amplitude  perturbations  ln(A/Ao)  at  10  seconds  for  finite  difference  calculation  (left)  and  Born 
approximation  (right).  The  star  is  at  the  source  location. 

Although  there  are  similarities  between  the  finite  difference  solution  and  Bom  approximation, 
they  clearly  do  not  agree  in  detail,  and  so  correcting  the  amplitude  at  a  particular  point  may  not 
give  a  more  accurate  result  than  an  uncorrected  amplitude,  particularly  at  higher  frequencies.  For 
example,  between  600  and  800  km  north  near  the  left  edge  of  the  grid,  the  Bom  approximation 
predicts  strong  signal  enhancement,  while  the  finite  difference  calculation  shows  amplitude 
reduction.  Figure  3.17  shows  snapshots  of  the  vertical  motion  at  6  times,  showing  clearly  how 
the  wavefront  bends  irregularly  and  interferes  with  itself  as  it  propagates.  The  dominant  period 
band  in  the  snapshots  is  10-20  seconds. 


29 


Norlti  [km)  North  [km)  North  [km) 


Heterogeneous  model.  attimeSOsec 


-23 


X  10 


200  400  600  800  1000  1200  1400 

East  (km) 

Heterogeneous  model.  atlimeSOOsec  ^  1o''^ 


200  400  600  800  1000  1200  1400 

East  (km) 


Heterogeneousmodel.attime1205ec  jjio' 


200  400  600  800  1000  1200  1400 

East  (km) 

Heterogeneous  model.  attime2405ec  y  iq' 


200  400  600  800  1000  1200  1400 


East  (km) 


Figure  3.17.  Snapshots  of  vertical  velocity  propagating  from  Lop  Nor  to  KNET  for  3D  finite  difference 

calculation.  Figures  show  snapshots  at  multiples  of  60  seconds  from  60  to  360  seconds.  Amplitude 
scale  range  is  proportional  to  1/t^^^,  which  approximately  corrects  for  surface  wave  geometric 
spreading. 
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4.0  Inversion  for  Earth  Structure  and  Attenuation 


4.1  The  inversion  procedure  for  a  3D  earth  model 


In  our  previous  projects,  we  inverted  a  large  volume  of  dispersion  data  for  global  earth  structure. 
Global  earth  structure  refers  to  a  set  of  vertically  layered  earth  models  defined  for  each  cell  of  a 
one-degree  by  one-degree  grid  of  the  earth.  This  procedure  is  summarized  here,  and  in  the 
following  sections  we  show  how  it  is  modified  to  include  scattering  and  diffraction  and  modified 
to  invert  attenuation  data  for  global  Q  structure.  The  relationship  between  dispersion  and  the 
shear  wave  velocities  of  the  layers  in  the  earth  model  is  non-linear,  so  the  shear  velocities  are 
estimated  by  an  iterative  least  squares  inversion  procedure.  At  each  step  a  system  of  equations  is 
formed,  augmented  by  additional  equations  of  constraint,  and  then  solved  by  the  LSQR 
algorithm  (Nolet,  1987).  The  equations  solved  are 
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Ax  = 

-sHx 
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(4.1) 


where  Ax  is  the  vector  of  adjustments  to  the  shear  wave  slownesses  of  layers  in  each  of  the  577 
model  types.  Ad  is  the  vector  of  slowness  differences  between  predicted  and  observed  dispersion 
measurements,  s  is  the  vector  of  residuals  that  remain  after  inversion  (the  inversion  minimizes 
|8|).  X  is  the  vector  of  slownesses  estimated  in  the  previous  iteration.  The  elements  of  the  matrix 
A  consist  of  partial  derivatives  of  dispersion  predictions  with  respect  to  shear  wave  slownesses  in 
each  layer.  H  is  a  difference  operator  that  applies  to  vertically  neighboring  layers  and  has  the 
effect  of  constraining  the  vertical  smoothness  of  velocity  profde.  H  applies  to  layers  in  the  crust 
and  upper  mantle,  but  has  explicit  discontinuities  at  the  crust/mantle  boundary  and  at  the  base  of 
surface  sediments,  s  is  the  weighting  of  the  smoothness  constraint  and  can  be  a  diagonal  matrix 
(for  variably  weighted  smoothing)  or  a  scalar  (constant  smoothing).  I  is  the  identity  matrix  and  X 
weights  the  damping  which  constrains  the  norm  of  the  difference  between  final  slownesses  and 
constraining  model  slownesses  Xc  (in  this  case  a  variant  of  the  Crust  2.0  values).  X  can  be  a  scalar 
for  constant  damping,  or  a  diagonal  matrix  for  variable  damping. 


4.2  Correction  for  scattering  and  diffraction  due  to  a  realistic  heterogeneous  earth  model 

The  Born  approximation  techniques  discussed  in  section  3  provide  a  straightforward,  but 
approximate,  way  to  incorporate  scattering  and  diffraction  into  the  inversion  procedure.  As 
described  above,  the  matrix  A  in  equation  4. 1  is  calculated  using  a  path  integral  to  calculate  the 
phase  velocity,  with  each  element  of  the  matrix  corresponding  to  a  piece  of  the  path  weighted 
according  to  the  fraction  of  the  path  crossing  a  grid  block  and  the  sensitivity  of  the  observable  to 
the  model  velocity.  This  can  be  replaced  by  integration  over  the  Fresnel  zone  area,  which 
changes  the  weighting  of  each  element  and  increases  the  number  of  elements  corresponding  to 
each  ray.  The  matrix  requires  more  time  to  calculate,  but  the  inversion  procedure  is  the  same  as 
in  the  ray-based  tomographic  inversion.  We  use  this  approach  for  inversion  of  dispersion  data  for 
earth  structure.  For  amplitudes,  we  use  the  earth  model  determined  by  inversion  of  dispersion 
data  to  calculation  amplitude  corrections  using  the  Bom  approximation  and  then  apply  these  to 
the  data  prior  to  inverting  for  Q  stmcture,  which  is  then  done  using  path  integrals. 
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4.3  Application  of  corrections  to  inversion  for  earth  structure 

We  ran  our  earth  structure  inversion  using  Bom  corrections  for  scattering  and  diffraction  by 
incorporating  the  finite  frequency  sensitivity  kernel  (equation  3.3)  into  the  inversion  code.  We 
then  reran  the  entire  global  tomographic  inversion  of  approximately  one  million  dispersion 
measurements  (Stevens  et  ah,  2005)  to  generate  a  new  set  of  earth  stractures.  We  found  that  the 
changes  from  the  previous  inversion  are  small  in  most  areas,  and  there  is  no  significant 
improvement  in  data  fit,  so  it  is  not  clear  that  the  results  represent  an  improvement  over  the 
inversion  using  great  circle  paths.  Because  the  differences  were  so  small,  we  used  the  existing 
stmctures  as  the  base  models  for  Q  inversion. 


4.4  Inversion  of  attenuation  data  for  Q  structure 

Inversion  of  attenuation  data  for  Q  stmcture  can  be  accomplished  using  equation  4.2,  which  has 
the  same  form  as  equation  4.1  above: 


( 

sH 

Ax  = 

-sHx 

^XI) 

X{xc-x) 

(4.2) 


with  the  following  changes: 

1.  The  data  are  attenuation  residuals  instead  of  dispersion  residuals.  Attenuation 
estimates  are  derived  from  an  existing  Q  model,  and  the  differences  between  those 
and  the  observations  are  the  data  used  in  the  inversion.  Amplitude  measurements  may 
be  corrected  for  the  effects  of  heterogeneous  stmcture. 

2.  The  matrix  A  is  derived  from  derivatives  of  the  attenuation  coefficients  with  respect 
to  model  Q  in  each  layer  for  the  path-averaged  inversion,  and  includes  Bom 
scattering  sensitivity  functions  for  the  area  integrals. 

3.  The  starting  model  and  constraining  model  are  the  same,  and  are  derived  from  PREM 
for  depths  greater  than  100  km,  and  Swanger’s  Law  (Q=100  [3,  with  [3  in  km/s)  for 
shallower  depths.  The  values  derived  from  PREM  are  Q=18  |3  for  depths  between  100 
km  and  220  km,  Q=30  [3  at  greater  depths.  There  are  discontinuities  at  100  km  and 
220  km  and  a  smoothness  criterion  is  applied  to  layers  above  and  below  100  km.  The 
inversion  is  performed  for  layers  shallower  than  220  km.  Q  is  fixed  to  1 8  p  at  220  km 
depth  and  to  30  p  below  this  depth. 

4.  The  model  vector  consists  of  p/Q  for  each  layer  that  is  free  to  change  in  each 
stmcture,  and  optionally  can  include  the  change  in  moment  of  each  event.  That  is,  for 
attenuation  residuals  that  were  derived  using  spectra  with  a  fixed  model  moment,  the 
moment  for  each  event  can  be  corrected  as  part  of  the  inversion. 

Unlike  inversion  for  shear  velocity,  the  inversion  for  Q  is  linear,  so  only  a  single  iteration  is 
necessary,  although  multiple  inversions  are  done  with  different  damping  and  smoothing 
parameters  to  generate  realistic  models.  Note  that  equations  3.1  and  3.2  give  the  equation  for  the 
predicted  surface  wave  spectral  amplitudes  given  a  source  mechanism  for  the  event.  Equation  3.4 
can  optionally  be  used  to  correct  for  stmctural  heterogeneity  using  the  Bom  approximation. 
Collecting  the  distance  independent  terms,  equations  3.1  and  3.2  can  be  written  in  the  form: 
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(4.3) 


I  .  ^1  S{co)A{a))ex^[-Yp{a))r] 

\uXo),r)\  =  Mo -  .  ,  ,  ' - 

^a^sin(r/aj 

where  Mg  is  the  source  moment  and  S(co)  the  source  function  (important  for  large  events),  A((o) 
is  a  frequency  dependent  function  that  depends  on  source  and  receiver  structure  and  focal 
mechanism.  We  assume  that  A((o)  can  be  predicted  well  enough  from  the  background  earth 
structure  and  the  source  mechanism  which  is  either  a  point  explosion  or  CMT  moment  tensor. 
While  the  inversion  program  has  the  capability  to  allow  an  amplitude  scale  factor  for  each  path, 
which  would  allow  to  variations  in  explosion  amplitude  due  to  tectonic  release,  and  variations  in 
earthquake  amplitude  due  to  errors  in  source  mechanism,  in  the  inversions  that  follow  we  only 
allowed  the  moment  to  vary,  which  is  a  constant  factor  for  all  paths  for  a  single  event.  In 
equation  4.3,  Mg  and  jp  are  derived  from  a  starting  source  mechanism  and  background  earth 
model,  and  allowed  to  vary  in  the  inversion  while  the  other  factors  are  held  fixed.  S(co)  is  derived 
for  a  triangular  function  with  rise  time  (half-duration)  T.  Since  this  is  approximate,  points  where 
fT  >  0.5,  which  corresponds  to  an  amplitude  reduction  of  0.4,  are  zero  weighted.  The 
observed  data  can  then  be  written: 


I  ^1  S(<D)A(co)cxp[-r;(<o)r] 

\u;  (®,  r)  =  M\ -  .  r  I  \ - 

^a^sm(r/aj 

And  the  log  ratio  of  observed  to  predicted  spectra  has  the  form: 

In  _  In  °  +  [y „ ico)rl  =  AlnM„  -l-rAy 


(4.4) 


(4.5) 


Ay  can  be  further  expanded  into  a  sum  over  Q  structure  in  each  structure  traversed  by  the  ray 
along  the  source  to  receiver  point  multiplied  by  the  fraction  of  ray  over  each  structure.  So  we  can 
rewrite  equation  4.5  as: 


1  ^  u;((o,r) 
r.,  u^ico,r) 


(4.6) 


The  subscript  k  refers  to  an  event,  i  refers  to  a  single  path  for  event  k,  1  refers  to  each  model  type 
traversed  by  the  ray  and  is  the  fraction  of  the  path  that  the  ray  spends  in  each  model.  The 

subscript  j  refers  to  layer  number  in  each  model.  L  is  the  total  number  of  model  types  traversed 
and  J  is  the  total  number  of  layers  allowed  to  change  in  each  model.  So  the  data  in  the  inversion 
is  the  left  hand  side  of  equation  4.6,  the  spectral  ratio  divided  by  the  distance  for  multiple 
frequencies,  and  the  inversion  is  performed  for  the  quantities  A  In  Mq  ,  the  change  in  moment  of 

each  event,  and  A(y^/g)^. ,  the  ratio  of  shear  velocity  to  Q  in  each  layer  of  each  model.  The 

function  Gij  gives  the  change  in  y  with  respect  to  change  iny^/gand  can  be  written  assuming  no 
bulk  attenuation  in  terms  of  partial  derivatives  of  phase  velocity  with  respect  to  material 
velocities  as: 


G  -  I  4  A- 

^Pij  3  a,j  da,j 


(4.7) 
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The  inversions  described  below  all  invert  data  of  the  form  described  above.  We  also  compare 
with  interstation  attenuation  estimates.  Using  equation  4.5  again  at  two  or  more  stations  at 
distances  rn,  we  get  a  set  of  equations  of  the  form: 


In 


u^{co,rJ 


=  A  hi  Mg  + 


(4.8) 


The  equations  define  a  line  with  slope  Ay  and  intercept  AlnMo  so  this  can  be  used  as  a  check  on 
the  inversion  by  finding  the  average  change  in  attenuation  over  all  or  a  subset  of  paths  for  each 
event,  and  the  change  in  moment. 


4.5  Data  used  in  the  Q  Inversion 

We  used  two  sets  of  data  in  our  Q  inversion.  The  first  data  set  was  provided  to  us  by  Anatoli 
Levshin  of  the  University  of  Colorado  and  is  described  in  Levshin  et  al.  (2007).  They  used 
essentially  the  same  procedure  described  above,  determining  gamma  from  equation  4.8,  except 
that  they  allowed  the  moment,  depth  and  fault  orientation  to  vary  in  order  to  determine  more 
realistic  attenuation  coefficients.  The  second  set  of  data  consisted  of  our  own  measurements  on  a 
different  set  of  data  also  covering  the  Eurasian  continent.  The  events  processed  are  listed  in 
Appendix  A  and  the  events  and  recording  stations  are  shown  in  Figure  4.1.  Figure  4.2  shows  a 
histogram  of  path  lengths  in  the  complete  data  set.  As  discussed  above,  we  determined 
attenuation  coefficients  using  CMT  moments,  but  then  allowed  the  moments  to  vary  as  part  of 
the  inversion  process.  There  is  quite  a  lot  of  scatter  in  both  data  sets,  and  we  rely  on  data 
redundancy  to  help  define  the  attenuation  model. 


Events 


Stations 


Figure  4.1  Maps  showing  events  processed  (left)  and  stations  recording  data  from  these  events  (right). 


Histogram  of  Path  Length  Distribution 


6000  SDOO 

Distance 


Figure  4.2.  Histogram  of  distances  of  all  rays  in  the  data  set.  Each  frequency  is  considered  distinct  in  this  plot. 
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Four  separate  inversions  were  performed  with  the  following  data  sets:  1)  SAIC  data,  2)  SAIC 
data  with  Bom  corrections  applied,  3)  CU  data,  and  4)  CU  and  SAIC  data  combined.  Figure  4.3 
shows  the  change  in  moment  determined  as  part  of  the  inversion.  Since  the  CU  data  was  moment 
corrected  prior  to  inversion,  the  change  in  moment  for  that  data  set  is  quite  small.  The  variation 
in  moment  for  the  SAIC  data  is  similar  to  that  described  by  Levshin  et  al.  (2007). 


Moment  Corrections  Derived  from  Q  Inversion 


Figure  4.3.  Change  in  moment  for  each  event  processed.  ^‘All  Data”  refers  to  inversion  of  the  SAIC  and  CU 
data  simultaneously  while  the  other  data  sets  were  inverted  independently.  Inversion  of  all  data 
together  gives  very  similar  results  to  inversion  of  individual  data  sets.  Application  of  Born 
corrections  also  makes  little  difference  to  the  moment  corrections. 

4.6  Q  Inversion  Results 

Inversion  results  were  evaluated  by  examining  the  models  and  the  data  fit.  The  models  were 
initially  inverted  with  the  same  damping  and  smoothing  factors  for  all  models.  Then  the  post¬ 
inversion  models  were  examined  to  look  for  problems.  The  most  common  problem  is  negative  Q 
values  which  are  a  sign  of  underdamping  of  the  inversion.  The  damping  and  smoothing 
parameters  were  both  increased  by  the  same  amount  for  these  models  and  the  inversion  rerun 
until  models  were  all  physically  reasonable.  Figure  4.4  shows  an  example  of  the  initial  model 
and  final  models  for  each  data  set. 


Figure  4.4  Comparison  of  inversion  results  with  base  model  for  one  central  Asian  structure.  In  this  case,  all 
data  sets  show  lower  Q,  with  the  Born  adjusted  data  set  showing  a  larger  change. 
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The  second  check  on  the  inversion  was  examination  of  data  fit.  This  was  done  both  for 
individual  models  and  for  averages  over  all  data.  Because  of  the  large  variability  in  the  data,  we 
typically  get  good  data  fits  for  some  paths  for  each  event,  and  some  that  are  not  in  good 
agreement  with  the  data,  however  for  each  event  we  generally  have  good  data  fits  for  most  paths. 
Figure  4.5  shows  examples  for  a  few  paths  for  one  event.  Figure  4.6-4.8  show  examples  for 
averages  along  many  paths. 


Figure  4.5.  Data  fits  for  four  paths  of  different  ranges  for  the  same  event.  The  red  line  is  the  starting  model 
and  the  green  line  the  final  gamma  model.  The  blue  line  corresponds  to  the  data  fit  and  is  the  green 
line  shifted  by  an  amount  corresponding  to  the  moment  adjustment.  This  adjustment  is  larger  at 
closer  distances.  The  red  marks  are  the  data  points. 

Figure  4.6  shows  the  average  of  the  data  over  all  paths  longer  than  5000  km,  the  inversion  results 
for  the  same  paths,  and  the  starting  model  for  the  same  paths.  Both  the  data  and  the  inversion 
results  show  substantially  higher  attenuation  than  in  the  starting  model.  The  average  inversion 
results  are  a  very  good  fit  to  the  average  data.  The  results  suggest  that  a  background  model  with 
Q=70p  in  the  crust  and  upper  mantle  would  be  more  consistent  with  the  data  than  our  starting 
model  with  Q=100p. 
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Gamma,  path  range  5DDD  10000 


Frequency  (Hz) 


Figure  4.6.  Average  gamma  for  all  paths  after  inversion  (blue),  in  original  data  set  (red)  and  starting  model 
(green).  Dashed  lines  are  ±1  standard  deviation.  For  paths  between  5000  and  10,000  km. 

Figure  4.7  shows  the  average  over  paths  shorter  than  5000  km  of  the  data,  the  inversion  results 
for  the  same  paths,  and  the  starting  model  for  the  same  paths.  The  attenuation  is  higher  along 
these  shorter  paths,  but  the  scatter  in  the  data  is  also  much  larger,  as  indicated  by  the  wide 
standard  deviations  on  the  plot. 


Figure  4.7.  Average  gamma  for  all  paths  after  inversion  (blue),  in  original  data  set  (red)  and  starting  model 
(green).  Dashed  lines  are  ±1  standard  deviation.  For  paths  between  1000  and  5,000  km. 

Figure  4.8  shows  the  average  over  all  paths  longer  than  1000  km  and  for  all  paths  longer  than 
5000  km  for  the  inversion  results  and  the  starting  model  for  the  same  paths.  Now  we  also 
compare  with  the  observed  amplitude  decay  along  all  paths  for  each  event.  The  attenuation 
determined  from  amplitude  decay  is  slightly  smaller  than  determined  from  the  inversion  results, 
but  well  within  the  scatter  in  the  data,  and  both  are  substantially  higher  than  the  starting  model, 
particularly  at  higher  frequencies. 
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Gamma  Compared  with  Residuals 


Frequency  (Hz) 


Frequency  (Hz) 


Figure  4.8.  Average  gamma  for  all  paths  after  inversion  (blue),  in  starting  model  (green),  and  as  determined 
by  ampiitude  decay  aiong  ail  paths  for  each  event  (red).  Dashed  lines  are  ±1  standard  deviation. 
Average  gamma  and  starting  model  are  for  paths  in  the  data  set  ionger  than  1000  km  (top)  and 
ionger  than  5000  km  (bottom). 

Figures  4.9-4.14  show  the  inversion  results  at  frequencies  of  0.05,  0.067,  0.08,  0.1,  0.125  and 
0.15  Hz.  Each  plot  shows  a  map  of  inversion  results  in  Eurasia  and  a  histogram  of  attenuation 
values  for  each  frequency.  The  results  from  the  different  data  sets  differ  in  some  details,  but  are 
similar  in  most  respects.  All  show  a  band  of  high  attenuation  stretching  across  Asia  from  the 
Middle  East  through  Southeast  Asia. 
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Gamma  (km'^}  0.05  Hz  -  All  Data 
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Gamma  (km^^)  0.06  Hz  -  SAIC  Born  Corrected  Data 
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Gamma  (km'^)  0  05  Hz  -  CU  Data 


Figure  4.9.  Attenuation  inversion  results  at  20  seconds  -  left:  gamma  map,  right:  gamma  histogram. 

Inversion  data  sets  from  top  to  bottom:  All  data,  SAIC  data,  SAIC  data  Born  Corrected,  CU. 
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Gamma  0.067  Hz  -  All  Data 
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Gamma  0.067  Hz  -  SAIC  Data 


Gamma  0.067  Hz  -  SAIC  Born  Corrected  Data 
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Figure  4.10  Attenuation  inversion  results  at  15  seconds  -  left:  gamma  map,  right:  gamma  histogram. 
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Gemma  0.06  Hz  -  All  Dp;g 
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Figure  4.11.  Attenuation  inversion  results  at  12.5  seconds  -  left:  gamma  map,  right:  gamma  histogram. 
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Gemma  0.10  Hz  -  All  Dplg 
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Gamma  (km'^)  0.10  Hz  -  All  Data 


Gamma  (km'^)  0.10  Hz  -  SAIC  Data 
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Gamma  (km'^)  O.'IO  Hz  -  SAIC  Born  Corrected  Data 
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Figure  4.12.  Attenuation  inversion  results  at  10  seconds  -  left:  gamma  map,  right:  gamma  histogram. 
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Gamma  (km"’)  0.10  Hz  -  All 
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Gamma  (km'^)  0.125  Hz  -  SAIC  Born  Correcl^d  Data 
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Gamma  (km  ')  0.125  Hz  -  CU  Data 


Figure  4.13.  Attenuation  inversion  results  at  8  seconds  -  left:  gamma  map,  right:  gamma  histogram. 
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Gemma  0.1^  Hz  -  All  Dp|g 
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Gamma  [km'^)  0.15  Hz  -  SAIC  Born  Corrected  Data 


-3 


x  10 


60  E 


Gamma  (km"")  0,15  Hz  ■  CU  Data 


-3 

X  10 


Figure  4.14.  Attenuation  inversion  results  at  6.6  seconds  -  left:  gamma  map,  right:  gamma  histogram. 
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5.0  Conclusions  and  Recommendations 

We  have  performed  a  detailed  study  of  high  frequency  surface  wave  amplitudes  and  attenuation 
in  Eurasia.  Following  are  the  major  sections  of  this  study  with  conclusions  and 
recommendations : 

1.  Regionalized  magnitudes  -  The  Russell  Butterworth  filtered  magnitude  and  the  path 
corrected  time  and  frequency  domain  magnitudes  have  a  similar  purpose,  specifically  to 
allow  surface  waves  to  be  measured  at  regional  and  local  distances  and  at  higher 
frequencies  than  20  seconds.  These  magnitudes  give  similar  results  when  applied  to  data 
and  are  consistent  in  value  with  traditional  20  second  magnitudes.  The  path  corrected 
magnitudes  have  the  advantage  that  they  can  be  regionalized  to  take  into  account 
differences  in  earth  structure  and  attenuation,  while  the  Butterworth  filtered  magnitude 
uses  a  good  representative  average  value  for  these  quantities.  An  issue  with  all  of  the 
magnitudes  is  how  to  determine  which  frequency(ies)  to  measure.  The  path  corrected 
spectral  magnitude,  for  example,  performs  a  robust  average  over  all  frequencies,  while 
the  Bonner  et  al.  implementation  of  the  Butterworth  filtered  magnitude  uses  the 
maximum  value.  Using  the  maximum  value  may  give  better  discrimination  but  at  the  cost 
of  more  variability  in  the  magnitude.  More  research  is  needed  to  determine  the  optimum 
procedure. 

2.  North  Korean  surface  wave  magnitude  -  The  North  Korean  nuclear  test  had  a 
surprisingly  large  surface  wave  magnitude,  nearly  a  magnitude  unit  higher  than  would  be 
expected  based  on  larger  events  assuming  a  Mgdog  yield  slope  of  one.  Some  of  this 
difference  can  be  explained  by  the  absence  of  tectonic  release  and  high  velocity  source 
medium  for  this  event.  However,  the  results  also  suggest  that  the  Mjdog  yield  slope  may 
be  less  than  one,  which  has  implications  for  discrimination  of  small  events  as  well  as 
yield  estimation.  More  research  is  needed  to  analyze  surface  waves  from  small  events  of 
known  yield  in  high  velocity  media  to  see  if  the  North  Korean  event  is  an  anomaly  or  the 
norm. 

3.  Bom  approximation  -  We  modeled  the  effect  of  heterogeneous  stracture  using  the  Bom 
approximation.  We  performed  both  data  analysis  and  large  3  dimensional  finite 
difference  calculations  to  assess  the  performance  of  the  Bom  approximation.  Our  goal 
was  to  use  the  Bom  approximation  to  correct  for  scattering  and  diffraction  caused  by 
heterogeneous  stmcture  prior  to  performing  Q  inversions.  However,  the  stmctural 
complexity  appears  to  exceed  the  limits  of  the  Born  approximation,  and  application  of  the 
Bom  corrections  do  not  improve  inversion  results  at  these  high  frequencies.  Amplitude 
variations  due  to  stmcture  at  these  high  frequencies  are  quite  large,  so  a  better  way  to 
correct  for  them  is  needed. 

4.  Q  inversions  -  We  inverted  surface  wave  amplitude  data  for  attenuation  and  corrections 
to  source  moment  for  data  from  about  300  Eurasian  earthquakes.  We  used  two  data  sets: 
one  our  own  measurements  and  one  set  from  Anatoli  Levshin  at  the  University  of 
Colorado.  The  data  sets  are  fairly  consistent  and  both  indicate  higher  attenuation  than  our 
initial  background  model.  There  is  a  large  amount  of  scatter  in  both  data  sets,  likely  the 
result  of  stmctural  variations  and  interference  as  discussed  above.  Nevertheless,  there  is 
enough  redundancy  in  the  data  that  we  were  able  to  perform  Q  inversions  for  the 
Eurasian  continent.  The  inversions  could  be  significantly  improved  in  two  ways:  1)  by 
implementing  a  better  correction  for  variations  due  to  earth  stmcture;  and  2)  by 
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increasing  the  size  of  the  data  set  and  measuring  attenuation  from  more  events.  The  latter 
is  straightforward,  but  the  former  is  still  troublesome.  Ray  tracing,  for  example,  is 
another  approach,  but  the  Bom  approximation  reduces  to  ray  tracing  in  the  high 
frequency  limit,  and  it  is  subject  to  similar  problems.  A  multiple  scattering  Bom 
approximation  as  suggested  by  Friederieh  and  Maupin,  would  be  another  possibility,  but 
it  is  not  known  whether  it  would  lead  to  better  results.  Finite  difference  calculations  such 
as  we  did  in  this  study  are  another  possibility,  but  they  take  considerable  eomputational 
time  and  require  detailed  knowledge  of  earth  stmcture.  The  most  promising  approach  is 
likely  to  be  a  hybrid  observational/eomputational  method  that  uses  the  data  together  with 
a  numerical  or  semi-empirical  model  for  surface  wave  propagation  to  match  observed 
waveforms. 

5.  Data  and  model  distribution  -  The  attenuation  models  developed  during  this  project  have 
been  incorporated  into  the  global  earth  models  in  the  format  distributed  earlier  by  Stevens 
et  al.  (2005)  and  are  available  to  all  researchers  with  the  permission  of  AFRL. 
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6.0  Data  Deliverable 


This  report  is  accompanied  by  a  data  deliverable  consisting  of  the  following  files: 

LP_2008_May.tar.gz:  This  is  the  complete  set  of  global  earth,  dispersion  and  attenuation 
models.  The  contents  are  described  in  the  file  “README.models.”  It  also  contains  a  program 
compiled  under  Linux  to  calculate  the  dispersion  and  attenuation  between  any  two  points  at  an 
input  frequency.  The  attenuation  models  are  the  “All  Data”  models  described  above  that  were 
derived  using  both  SAIC  and  University  of  Colorado  attenuation  measurements. 

gamma_data:  These  are  all  of  the  attenuation  measurements  made  by  SAIC  during  this  project. 
The  format  of  the  file  is  given  in  the  file  “dataformaf ’. 

moments.txt:  This  file  gives  the  adjustments  to  CMT  moments  found  during  the  Q  inversion 
process  as  described  in  Section  4  of  this  report.  The  columns  correspond  to  event  number, 
natural  log  of  the  correction,  and  value  of  the  correction. 

Note  that  the  data  in  the  “gamma_data”  file  are  the  original  measurements  derived  using  CMT 
moments  and  so  they  should  be  adjusted  using  the  moment  corrections  by: 

Ay  =  -AlnMo/r  (6.1) 

This  means  that  if  the  inversion  showed  that  the  moment  is  larger  than  the  CMT  moment  then 
gamma  should  be  reduced  for  each  data  point,  and  if  the  moment  was  found  to  be  smaller,  then 
gamma  should  increase. 
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Appendix  A.  Events  processed  for  Q  inversion 

The  following  table  lists  the  events  that  were  processed  at  SAIC  for  Q  inversion.  Evid  is  a 
sequence  number,  Date/Time  the  date  and  time  of  the  event,  Lat,  Lon  and  depth  the  event 
hypocenter,  Moment  is  the  CMT  moment  and  Rtime  the  rise  time  of  a  triangular  source  function, 
also  from  the  CMT  catalog,  mb  and  Mg  are  body  and  surface  wave  magnitude  and  Nsta  is  the 
number  of  stations  processed  for  the  event. 


Table  Al.  List  of  events  processed  for  Q  inversion. 


EVID 

Date/Time 

Lat 

Lon 

Depth 

Moment 

Rtime 

mb 

Ms 

Nsta 

1 

1994/01/11  00:51:59 

25.23 

97.20 

33 

1.60E+18 

2.7 

5.9 

5.9 

14 

2 

1994/02/23  08:02:05 

30.85 

60.60 

10 

1.72E+18 

2.9 

6.0 

6.1 

15 

3 

1994/02/24  00:11:12 

30.78 

60.50 

13 

3.30E+18 

3.7 

6.0 

6.1 

12 

4 

1994/02/26  02:31:11 

30.90 

60.55 

12 

1.39E+18 

2.5 

5.8 

5.9 

16 

5 

1994/03/01  03:49:01 

29.10 

52.62 

17 

1.37E+18 

2.8 

5.8 

6.0 

13 

6 

1994/04/06  07:03:27 

26.19 

96.87 

33 

7.15E+17 

2.0 

5.6 

5.6 

19 

7 

1994/04/13  04:00:51 

22.78 

123.63 

36 

5.76E+17 

2.1 

5.7 

5.6 

5 

8 

1994/05/01  12:00:37 

36.90 

67.16 

26 

1.65E+18 

2.6 

5.9 

6.3 

14 

9 

1994/05/23  05:36:03 

24.17 

122.54 

34 

1.89E+18 

2.9 

5.7 

6.0 

10 

10 

1994/05/23  15:16:58 

24.07 

122.56 

33 

8.13E+17 

2.2 

5.9 

5.7 

11 

11 

1994/05/24  04:00:46 

23.96 

122.45 

47 

6.60E+18 

4.2 

6.0 

6.6 

9 

12 

1994/05/29  14:11:51 

20.56 

94.16 

42 

6.50E+18 

4.4 

6.2 

6.2 

6 

13 

1994/06/05  01:09:31 

24.51 

121.91 

16 

3.80E+18 

3.8 

6.0 

6.5 

9 

14 

1994/06/20  09:09:04 

28.97 

52.61 

17 

8.05E+17 

2.1 

5.9 

5.7 

7 

15 

1994/06/29  18:22:36 

32.57 

93.67 

33 

7.72E+17 

2.2 

5.8 

5.5 

14 

16 

1994/08/19  21:02:45 

17.97 

96.42 

12 

4.81E+17 

1.8 

5.5 

5.6 

9 

17 

1994/08/21  15:56:01 

56.76 

117.90 

33 

1.25E+18 

2.6 

5.7 

5.8 

19 

18 

1994/09/16  06:20:18 

22.55 

118.74 

12 

1.25E+19 

5.4 

6.5 

6.7 

12 

19 

1994/11/21  08:16:36 

25.49 

96.70 

33 

9.25E+17 

2.4 

5.6 

5.9 

11 

20 

1995/02/23  21:03:02 

35.05 

32.28 

15 

8.06E+17 

2.0 

5.8 

5.7 

19 

21 

1995/02/23  05:19:02 

24.14 

121.61 

44 

2.45E+18 

3.5 

5.8 

6.2 

11 

22 

1995/05/13  08:47:12 

40.15 

21.70 

13 

7.64E+18 

4.3 

6.2 

6.5 

18 

23 

1995/06/15  00:15:48 

38.40 

22.28 

14 

6.01E+18 

4.3 

6.0 

6.5 

12 

24 

1995/06/25  06:59:05 

24.60 

121.71 

47 

1.02E+18 

2.4 

5.8 

5.7 

14 

25 

1995/06/29  23:02:31 

51.96 

103.10 

33 

5.20E+17 

1.8 

5.6 

5.5 

30 

26 

1995/07/09  20:31:31 

21.98 

99.16 

12 

7.53E+17 

2.1 

5.7 

5.9 

11 

27 

1995/07/11  21:46:39 

21.97 

99.20 

13 

1.91E+19 

6.6 

6.1 

7.2 

15 

28 

1995/10/01  15:57:16 

38.06 

30.13 

33 

4.72E+18 

4.0 

5.7 

6.1 

11 

29 

1995/10/23  22:46:50 

26.00 

102.23 

0 

2.18E+18 

3.3 

5.5 

0.0 

17 

30 

1995/11/13  08:43:14 

56.10 

114.50 

24 

5.50E+17 

1.8 

5.9 

5.6 

30 

31 

1995/11/22  04:15:11 

28.83 

34.80 

10 

7.21E+19 

11.0 

6.2 

7.3 

23 

32 

1995/12/05  18:49:33 

39.44 

40.15 

29.3 

5.54E+17 

2.0 

4.9 

0.0 

19 

33 

1996/02/03  11:14:19 

27.29 

100.28 

10 

9.94E+18 

5.2 

6.3 

6.5 

14 

34 

1996/03/05  14:52:28 

24.09 

122.22 

30 

3.59E+18 

3.9 

6.1 

6.4 

17 

35 

1996/03/05  17:32:10 

24.03 

122.24 

33 

7.33E+17 

2.3 

5.6 

5.6 

14 

36 

1996/03/19  15:00:26 

39.99 

76.70 

28 

3.60E+18 

3.9 

5.7 

6.0 

24 

37 

1996/03/29  03:28:56 

24.14 

122.20 

33 

5.20E+17 

1.8 

5.4 

5.5 

14 

38 

1996/05/03  03:32:47 

40.77 

109.66 

26 

1.07E+18 

2.6 

5.5 

6.0 

12 

39 

1996/06/22  16:47:12 

75.82 

134.62 

10 

4.95E+17 

1.6 

5.6 

5.5 

20 

40 

1996/07/20  00:00:41 

36.15 

27.10 

33 

2.38E+18 

3.0 

5.7 

6.2 

13 

41 

1996/08/14  02:59:41 

40.75 

35.34 

10 

3.54E+17 

1.6 

5.2 

5.5 

13 

42 

1996/09/05  23:42:06 

21.90 

121.50 

20 

1.91E+19 

6.4 

6.4 

6.6 

11 
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43 

1996/10/09  13:10:52 

34.56 

32.13 

33 

1.85E+19 

6.9 

6.4 

6.8 

12 

44 

1996/11/19  10:44:46 

35.35 

78.13 

33 

2.37E+19 

6.9 

6.1 

7.1 

26 

45 

1996/11/20  18:09:19 

39.60 

96.68 

33 

5.13E+17 

1.9 

5.0 

0.0 

2 

46 

1997/01/09  13:43:31 

41.03 

74.28 

22 

5.66E+17 

1.9 

5.7 

5.8 

40 

47 

1997/01/21  01:48:30 

39.47 

77.00 

33 

7.74E+17 

2.1 

5.3 

5.8 

25 

48 

1997/02/04  10:37:47 

37.66 

57.29 

10 

6.72E+18 

4.5 

5.9 

6.8 

29 

49 

1997/02/27  21:08:02 

29.98 

68.21 

33 

5.20E+19 

7.9 

6.3 

7.3 

9 

50 

1997/02/28  12:57:18 

38.08 

48.05 

10 

1.73E+18 

2.8 

5.5 

6.1 

18 

51 

1997/03/20  08:50:40 

30.14 

68.02 

33 

8.09E+17 

2.0 

5.5 

5.8 

34 

52 

1997/04/05  23:46:19 

39.51 

76.87 

33 

7.73E+17 

2.2 

5.4 

5.9 

28 

53 

1997/04/06  04:36:35 

39.54 

77.00 

33 

1.05E+18 

2.3 

5.6 

5.8 

35 

54 

1997/04/11  05:34:42 

39.53 

76.94 

15 

2.06E+18 

3.0 

5.8 

6.1 

35 

55 

1997/04/15  18:19:10 

39.63 

76.99 

23 

6.56E+17 

2.2 

5.4 

5.8 

28 

56 

1997/05/08  02:53:14 

24.89 

92.25 

35 

8.57E+17 

2.3 

5.6 

5.6 

18 

57 

1997/05/10  07:57:29 

33.83 

59.81 

10 

7.35E+19 

1.0 

6.4 

7.3 

24 

58 

1997/05/21  22:51:28 

23.08 

80.04 

36 

5.83E+17 

1.9 

6.0 

5.6 

28 

59 

1997/06/25  19:38:40 

33.94 

59.48 

10 

7.40E+17 

2.2 

5.5 

5.8 

29 

60 

1997/11/08  10:02:52 

35.07 

87.33 

33 

2.23E+20 

14.7 

6.2 

7.9 

43 

61 

1997/12/30  13:43:18 

25.38 

96.61 

33 

5.14E+17 

1.9 

5.4 

5.7 

18 

62 

1998/01/10  03:50:41 

41.08 

114.50 

30.3 

4.48E+17 

1.8 

5.8 

5.7 

29 

63 

1998/02/04  14:33:21 

37.08 

70.09 

33 

8.36E+17 

2.3 

5.6 

6.1 

38 

64 

1998/03/14  19:40:27 

30.15 

57.61 

9 

9.43E+18 

5.0 

5.9 

6.9 

26 

65 

1998/04/10  15:00:53 

32.46 

59.98 

33 

5.01E+17 

1.9 

5.3 

5.7 

38 

66 

1998/05/03  23:30:21 

22.31 

125.31 

33 

1.83E+20 

14.1 

6.4 

7.3 

14 

67 

1998/05/30  06:22:29 

37.11 

70.11 

33 

7.89E+18 

5.0 

5.9 

6.9 

34 

68 

1998/06/27  13:55:52 

36.88 

35.31 

33 

2.96E+18 

3.5 

5.8 

6.2 

34 

69 

1998/07/20  01:05:58 

30.13 

88.17 

33.2 

4.77E+17 

1.8 

5.4 

5.4 

62 

70 

1998/07/24  18:44:04 

21.25 

122.02 

33 

1.73E+18 

2.9 

5.6 

5.9 

23 

71 

1998/08/25  07:41:40 

30.08 

88.11 

33 

6.81E+17 

1.9 

5.3 

5.5 

51 

72 

1998/08/27  09:03:36 

39.66 

77.34 

33 

3.89E+18 

3.9 

5.6 

6.4 

42 

73 

1998/11/19  15:39:19 

22.61 

125.78 

10 

3.27E+18 

3.6 

5.8 

6.0 

15 

74 

1999/02/11  14:08:51 

34.26 

69.36 

33 

1.27E+18 

2.6 

5.4 

5.8 

41 

75 

1999/02/22  13:49:00 

24.12 

122.65 

42.7 

8.05E+17 

2.3 

5.4 

5.6 

27 

76 

1999/02/25  18:58:29 

51.60 

104.86 

10 

8.91E+17 

2.2 

5.9 

5.5 

30 

77 

1999/03/04  05:38:26 

28.34 

57.19 

33 

1.01E+19 

5.5 

6.2 

6.5 

38 

78 

1999/03/21  16:16:02 

55.90 

110.21 

10 

8.50E+17 

2.1 

5.5 

5.7 

39 

79 

1999/03/28  19:05:11 

30.51 

79.40 

15 

7.77E+18 

4.7 

6.4 

6.6 

52 

80 

1999/05/06  23:00:53 

29.50 

51.88 

33 

2.47E+18 

3.5 

5.9 

6.3 

24 

81 

1999/08/17  00:01:39 

40.75 

29.86 

17 

2.88E+20 

20.7 

6.3 

7.8 

18 

82 

1999/09/07  11:56:49 

38.12 

23.61 

10 

1.14E+18 

2.6 

5.6 

5.8 

12 

83 

1999/09/13  11:55:28 

40.71 

30.05 

13 

5.96E+17 

2.0 

5.8 

5.8 

42 

84 

1999/09/20  17:47:18 

23.77 

120.98 

33 

3.38E+20 

19.9 

6.5 

7.7 

40 

85 

1999/09/20  21:46:42 

23.39 

120.96 

33 

4.83E+18 

4.0 

5.8 

6.5 

37 

86 

1999/09/22  00:14:39 

23.73 

121.17 

26 

5.03E+18 

4.0 

6.2 

6.4 

48 

87 

1999/09/22  00:49:42 

23.64 

121.14 

33 

6.31E+17 

2.0 

5.9 

5.9 

38 

88 

1999/09/25  23:52:48 

23.74 

121.16 

17 

6.01E+18 

4.2 

6.2 

6.4 

38 

89 

1999/10/22  02:18:58 

23.45 

120.51 

33 

6.95E+17 

2.1 

5.7 

5.6 

39 

90 

1999/11/01  17:53:00 

23.38 

121.52 

33 

3.29E+18 

3.6 

6.1 

6.1 

49 

91 

1999/11/12  16:57:19 

40.76 

31.16 

10 

6.65E+19 

10.5 

6.3 

7.5 

17 

92 

1999/12/03  17:06:54 

40.36 

42.35 

19.3 

3.97E+17 

1.6 

5.3 

5.5 

24 

93 

2000/01/14  23:37:07 

25.61 

101.06 

33 

8.33E+17 

2.3 

5.4 

5.9 

54 

94 

2000/01/28  14:21:07 

43.05 

146.84 

61.1 

1.98E+19 

6.4 

6.7 

6.6 

28 

51 


95 

2000/06/06  02:41:49 

40.69 

32.99 

10 

1.11E+18 

2.7 

5.5 

6.1 

31 

96 

2000/06/07  21:46:55 

26.86 

97.24 

33 

3.74E+18 

3.7 

6.3 

6.5 

50 

97 

2000/06/10  18:23:29 

23.84 

121.23 

33 

5.35E+18 

4.3 

6.2 

6.2 

32 

98 

2000/08/04  21:13:02 

48.79 

142.25 

10 

1.92E+19 

6.4 

6.3 

7.1 

31 

99 

2000/08/22  16:55:13 

38.12 

57.38 

10 

3.62E+17 

1.7 

5.2 

5.8 

56 

100 

2000/09/10  08:54:46 

24.01 

121.53 

34.8 

5.83E+17 

2.0 

5.6 

5.6 

45 

101 

2000/09/12  00:27:58 

35.39 

99.34 

10 

1.76E+18 

3.1 

5.7 

6.3 

71 

102 

2000/10/09  02:30:00 

10.00 

92.95 

33 

3.57E+17 

2.0 

5.3 

5.6 

32 

103 

2000/12/06  17:11:06 

39.57 

54.80 

30 

3.90E+19 

8.8 

6.7 

7.5 

65 

104 

2000/12/15  16:44:47 

38.46 

31.35 

10 

1.21E+18 

2.7 

5.1 

5.8 

28 

105 

2001/01/03  14:47:49 

43.93 

147.81 

33 

7.66E+17 

2.1 

5.9 

5.2 

20 

106 

2001/01/26  03:16:40 

23.42 

70.23 

16 

3.43E+20 

24.1 

6.9 

8.0 

53 

107 

2001/01/28  01:02:10 

23.51 

70.52 

10 

5.22E+17 

2.0 

5.9 

5.5 

47 

108 

2001/03/05  15:50:07 

34.37 

86.90 

33 

8.64E+17 

2.1 

5.4 

5.8 

60 

109 

2001/03/15  01:22:43 

8.66 

94.01 

33 

1.04E+18 

2.3 

5.6 

5.9 

32 

110 

2001/03/24  06:27:53 

34.08 

132.53 

50 

1.97E+19 

6.2 

6.4 

6.5 

33 

111 

2001/06/10  01:52:08 

39.84 

53.89 

34.1 

1.49E+17 

1.2 

5.5 

5.2 

29 

112 

2001/06/14  02:35:25 

24.51 

122.03 

32.1 

7.80E+17 

2.2 

5.7 

5.6 

30 

113 

2001/06/25  13:28:46 

37.20 

36.17 

5 

1.68E+17 

1.1 

5.2 

4.9 

2 

114 

2001/07/26  00:21:36 

39.06 

24.34 

10 

5.61E+18 

4.0 

6.0 

6.6 

25 

115 

2001/11/14  09:26:10 

35.95 

90.54 

10 

5.90E+20 

25.8 

6.1 

8.0 

32 

116 

2001/12/18  04:02:58 

23.95 

122.73 

14 

2.08E+19 

7.1 

6.3 

7.3 

47 

117 

2002/02/03  07:11:28 

38.57 

31.27 

5 

6.00E+18 

4.4 

5.7 

6.4 

45 

118 

2002/02/03  09:26:43 

38.63 

30.90 

10 

6.11E+17 

1.7 

5.7 

5.6 

32 

119 

2002/02/12  03:27:25 

23.70 

121.57 

54.8 

3.79E+17 

1.8 

5.8 

5.4 

29 

120 

2002/02/17  13:03:52 

28.11 

51.76 

33 

1.16E+17 

1.1 

5.6 

5.0 

19 

121 

2002/03/25  14:56:33 

35.97 

69.17 

8 

1.62E+18 

3.4 

5.9 

6.2 

23 

122 

2002/03/27  08:52:52 

35.92 

69.28 

10 

2.83E+17 

1.5 

5.9 

5.4 

20 

123 

2002/03/31  06:52:50 

24.41 

122.21 

32.8 

5.45E+19 

10.1 

6.4 

7.4 

34 

124 

2002/04/12  04:00:23 

35.96 

69.42 

10 

7.24E+17 

2.4 

5.8 

5.9 

52 

125 

2002/04/17  08:47:22 

27.61 

56.76 

33 

9.11E+16 

1.0 

5.3 

4.9 

10 

126 

2002/04/24  10:51:50 

42.43 

21.51 

10 

4.50E+17 

1.9 

5.6 

5.6 

4 

127 

2002/05/15  03:46:05 

24.64 

121.92 

10 

1.91E+18 

3.0 

5.5 

6.2 

15 

128 

2002/05/28  16:45:17 

24.07 

122.26 

33 

1.49E+18 

2.5 

5.8 

5.9 

15 

129 

2002/06/04  14:36:05 

30.54 

81.44 

33 

2.99E+17 

1.6 

5.4 

5.3 

22 

130 

2002/06/22  02:58:21 

35.63 

49.05 

10 

6.97E+18 

4.6 

6.2 

6.4 

26 

131 

2002/07/1 1  07:36:26 

24.08 

122.29 

43.8 

6.52E+17 

2.1 

5.6 

5.6 

36 

132 

2002/07/13  20:06:27 

30.80 

69.98 

33 

4.78E+17 

1.9 

5.4 

5.7 

67 

133 

2002/08/08  1 1 :42:05 

30.99 

99.90 

33 

9.60E+16 

1.0 

5.4 

4.7 

15 

134 

2002/09/06  01:21:28 

38.38 

13.70 

5 

9.69E+17 

2.3 

5.8 

5.5 

13 

135 

2002/09/13  22:28:29 

13.04 

93.07 

21 

6.35E+18 

4.5 

6.2 

6.7 

37 

136 

2002/09/14  19:58:37 

13.06 

93.16 

33 

4.70E+17 

1.8 

5.7 

5.6 

33 

137 

2002/09/25  22:28:11 

32.09 

49.23 

10 

3.02E+17 

1.6 

5.5 

5.1 

8 

138 

2002/11/20  21:32:30 

35.41 

74.52 

33 

3.54E+18 

3.5 

5.7 

6.5 

49 

139 

2002/12/04  11:30:53 

19.40 

94.48 

53.5 

2.66E+17 

1.5 

5.6 

0.0 

11 

140 

2002/12/14  13:27:29 

39.73 

97.42 

22 

2.31E+17 

1.5 

5.6 

5.3 

11 

141 

2002/12/25  12:57:03 

39.70 

75.18 

10 

3.51E+17 

1.6 

5.5 

5.6 

24 

142 

2007/01/08  17:21:49 

39.80 

70.31 

14 

1.33E+18 

2.5 

5.9 

5.9 

33 

143 

2007/01/25  10:59:17 

22.56 

121.93 

36.2 

1.35E+18 

2.5 

5.6 

6.0 

21 

144 

2007/02/02  22:32:18 

37.71 

91.81 

10 

1.45E+17 

1.2 

5.3 

5.3 

29 

145 

2007/02/21  11:05:29 

38.43 

39.24 

10 

4.65E+17 

1.8 

5.7 

5.7 

20 

146 

2007/04/10  13:56:53 

13.00 

92.60 

30 

1.97E+17 

1.3 

5.5 

5.4 

16 

52 


147 

2007/05/05  08:51 :39 

34.25 

81.97 

9 

1.54E+18 

2.6 

5.8 

6.0 

4 

148 

2007/06/02  21:34:58 

23.03 

101.02 

10 

1.64E+18 

2.7 

6.2 

6.2 

19 

149 

2007/06/18  14:29:49 

34.48 

50.84 

10 

2.36E+17 

1.4 

5.5 

5.5 

29 

150 

2007/07/04  01:23:24 

55.50 

110.22 

10 

1.39E+17 

1.2 

5.3 

5.3 

12 

151 

2007/07/20  10:06:52 

42.93 

82.31 

11.8 

2.76E+17 

1.5 

5.5 

5.5 

18 

152 

2007/07/23  13:40:02 

23.64 

121.57 

40.6 

9.14E+16 

1.0 

5.6 

4.9 

25 

153 

2007/07/30  22:42:05 

19.31 

95.56 

14.2 

3.04E+17 

1.5 

6.0 

5.6 

8 

154 

2007/07/31  15:07:35 

27.31 

126.83 

10 

8.41E+17 

2.1 

5.5 

5.9 

8 

155 

2007/08/02  02:37:42 

47.12 

141.80 

5 

2.35E+18 

3.0 

5.3 

6.2 

29 

156 

2007/08/02  05:22:17 

46.71 

141.75 

10 

5.42E+17 

1.8 

5.6 

5.8 

23 

157 

2007/08/29  03:00:18 

21.73 

121.37 

24.5 

1.72E+17 

1.3 

5.7 

5.5 

10 

158 

2007/09/06  17:51:27 

24.33 

122.32 

62.9 

2.93E+18 

3.2 

6.5 

6.5 

16 

53 


